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Recursive Computation of the Repeated Integrals 
of the Error Function 


By Walter Gautschi 


1. This paper is concerned with a special technique, originated by J. C. P. 
Miller [1, p. xvii], of computing a solution f, of a second-order difference equation 


(1.1) Ynsi + GnYn + OnYn1 = O (n = 1,2,3,---) 


“» 


for n = 0(1)N, N large, in cases where (1.1) has a second solution, g, , which ul- 
timately grows much faster than f, . Straightforward use of (1.1) is then not ade- 
quate, since rounding errors will ‘“‘activate” the second solution g, , which in turn 
will eventually overshadow the desired solution f,. Miller’s device consists of 
applying (1.1) in backward direction, 


(1.2) Yn = —bn (nYn + Yn+i) , y eat a Ae Vv > N), 


starting with the initial values 
(1.3) Yi=a, y = 0, 


where a is any real number ~0. If v is taken sufficiently large the values so obtained 
turn out to be approximately proportional to f, in the range 0 < n < N. The factor 
of proportionality may then be determined, e.g., by comparing yo with fo . 

This technique was originally devised [1] for the computation of Bessel functions 
I,(x), and has since then been applied to various other Bessel functions [2], [5], 
[9], to Legendre functions [8], and to the repeated integrals of the error function* {6}, 

- -(t— og 
i* erie x = —. | ( t)" ay (n = 0,1,2, ---), 
VT z n! 
(1.4) 
2 eo 
Vr 


An analogous technique for first-order difference equations is described in [4, p. 25]. 

We shall present in Section 2 a detailed description of Miller’s procedure, 
paying special attention to the error term. In Section 3 we study the procedure as 
applied to the computation of the functions (1.4) and show that the process con- 
verges for any positive z, as vy — . In Sections 4-5 estimates will be developed of 
how large v must be taken to ensure any prescribed accuracy. 


—l 
t erfex = 


Received August 22, 1960. 
* In this notation i"(n = 0) should be interpreted as the nth power of the integral operator 


io) 
i-/ , so that 
z 


x 
® erfe z = erfe z, i* erfe x = i” erfc t dt (n = 1, 2, ---). 
Zz 


This notation for the repeated integrals of the error function, even though not entirely satis- 
factory, has become standard. 
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2. Consider the homogeneous second-order difference equation 


(2.1) Yat + AnYn + bnYn—1 =0 (n - i, 2, 3, ns +) 
and assume that 
(2.2) b. = 0 forall n 2 1. 


Let f, ve the (nontrivial) solution of (2.1) to be computed for n = 0(1)N. We 
assume 


(2.3) fo ¥ 0. 


Let there be another solution g, of (2.1), for which 





(2.4) gn 0 forall n 2 0, 
and 
| 
(2.5) lim | fn = 0. 
m-7eo | Gn | 


It follows readily that f, , gn are linearly independent. 

Now let yn” (n = 0,1, --- , v — 2; » > N) be the result of applying (2.1) in 
backward direction, starting with 
(2.6) yr=a yy =0 (a #0). 


These values, by (2.2), are well defined, and, as will presently be shown, yo” ¥ 0 
for v sufficiently large. Let us then define 


(2.7) tg = a Ya". 
We show that for any fixed n, 
(2.8) lim fa” = fa. 
Moreover, 
1 ioe f Jn 
(2.9) fo? = fy — BE. 
1 — Fe % 
gy fo 


It is sufficient to prove (2.9), since (2.8) then follows from (2.5). Let yn” be 
extended to all n > v by means of (2.1). Then for every fixed v the sequence 
{yn}(n = 0, 1, 2, ---) is a solution of (2.1), and therefore representable in the 
form 


— _ A” 4 + Bg, (n 


IV 


0). 
By (2.6), 
A”f,_, + B”g,_4 a, 
A” f, + Bg, = 0. 


(2.10) 





Th 


0). 
£0 


v) be 
ence 
1 the 


> 0). 
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Certainly, A” + 0, since otherwise, by (2.4), A” = B” = 0, which contradicts 
the first equation in (2.10). From the second equation, B” /A” = —f,/g,. Therefore 


4 ’ B” . ; 
Yn” - A” (1.4 5) =: A® (4, -£o,). 


If v is sufficiently large it follows because of (2.3) and (2.5) that yo” + 0. By (2.7), 
h A” (r. ~ £9.) 1 _ fr Gn 
(v) g 











fn - ’ -_ 2. gy , a 
A® (4 -& we) 1 _ Ff, % 
YJ» gy fo 
which proves (2.9). 
It is convenient to define 
(2.11) Pn = Jn 9% (n - 0, 1, 2, vee), 


Jn fo 
so that p, > 0 asn — ~, and 
() _ 1— (p»/pn) 
fn _— | A 2 . 
The relative error of f,'” is given by 
go a £. p ( 1 
2.12 = es & 
( Jn l—»p ) 


The approximations f,”” obviously do not depend on a, so that a can be chosen 
at will. If a high-speed computer is employed it is advisable to choose a small value 
for a to guard against “overflow” in the values of y,"”. 

3. Now let 





(3.1) Ie om i” erfe z, x 70 (n = 0, Re 2, vee) 
Then f, is a solution of 

x 1 7 
(3.2) Yat + n et an Yu = 0 (n _ 1, 2, 3, wi *) 


as is readily verified by writing 


n PS 2 1 -(¢—s)*" —{2 _@ -(¢—2z)*"* —t2 ) 
. eter = > (1 f a-nr”: “""st boa’ * 





and evaluating the first integral by parts. A second solution of (3.2) is given by 
(3.3) gn = (—1)"i”” erfe (—z) (n = 0,1, 2, ---). 
It is clear that the assumptions (2.2)—(2.4) are satisfied in this case. We shall now 


verify (2.5), i.e. 


: i" erfe x 
(3.4) — i* erfe (—2) =0 (x > 0). 
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This then will prove the convergence of the procedure in Section 2, as applied to 
(3.2). 

We recall that the repeated integrals of the error function are related to the 
parabolic cylinder functions D,(x) by [7, p. 76] 


—jz2 
é 


Q(n—1) 127 


It is furthermore known [8, p. 123] that 


D_,~(z) = Vre-¥ aller [1 +0 (+) (n — o, 2 bounded). 
9 (n+) /27, n Vn 
ro (3 + ') 


Therefore we obtain immediately for any fixed x, real or complex, 


miata at | 
ar 2 eric 7 = = an ae E +0 (2) (n— «), 
(3.5) op (2 ; Jn 
2 5) + 


i" erfe x = 


= D_»-a(arv/2 ). 
us 





Hence, 
i" erfe x 


————_—_——_ ~ ev 2nz (n a ) 
1” erfe (—2x) 


which proves (3.4). 
4. With f,, , gn defined by (3.1) and (3.3) we have for the quantities p, in (2.11) 


-n 


sh 
a erfex 
(n 


mG. a ey 
i”! erfe (—2)- ae? ) 


(4.1) pn = (—1)" 


It is shown in this section that for any fixed x > 0 the sequence {| p, |} is mono- 
tonically decreasing, i.e., 








| Pn+i | t” erfe x "ibe erfc (—2z) 
' ede _ —— <1 > 0). 
ae | Pn | “I erfex i" erfe (—2z) (n2 


Inequality (4.2) is obvious if n = 0 and, by (1.4), equivalent to 
/ (t — 2x)" e" dt | (s+2)"e ds 
ve / (t —2)"e" at | (s+2)"'e” ds>0 


if n > 0. By introducing new variables of integration, = u + 2,s = v — a, and 
writing the left-hand side as a double integral, one obtains 


(4.3) il] uty" (yn — wet -@-* du dv > 0, 
Q 


where Q denotes the first quadrant u = 0, v 2 0. Let Q: , Q2 denote the regions 


Q:: 0<u<y, Qe: 0<v<u. 





or 


hay 


whi 


Oak 
Oak 


hie 


no- 


> 0 


, and 


ons 
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Interchanging variables of integration gives 


n—1 n—1 —(u+z)2—(r- x)? 1 n—l —(v+z)2—(a )2 
If uv" (v—ue "~~ ”” dudv= -|/ uv (v—uje”~ "” dud. 
Qe Qi 


Therefore (4.3) is equivalent to 


_ = _ t+) 2—(y—z)2 a 2—(y4—r)2 
I u” _? Wy dd u)le (u+z) (v—z) _—* (v+z) (u—z) ] du dv > 0. 
Q 


Now, wu” 'v” “(v — u) > Oin Q, , and the same is true for the expression in brackets, 
since 
—(u+2)? —(v—2z) > —(v+2)* —(u—2z) for u<v. 
This proves (4.3), and thus (4.2). 
5. We are now in a position to estimate v such that for any given integer p, 


| + aa — fn) Te < 10” for n= 0, a ie N + B 


n 


Here, f,”’ denotes the approximations to f, = i" ' erfe x obtained by the procedure 
of Section 2. 
Since, by (2.12), 


(vy) 1 2 
(fn — fn)/fn < Py (1 + | Pn | ) _ 0( p, ), 
and since | p, | increases with n, by (4.2), it is sufficient to choose v such that 


(5.1) p»| (1+ | pwail) S 10°”. 


From (3.5) and (4.1) we have 


(5.2) | pasa | = cevene| + o( , )}. 
vn 


Assuming N large enough to neglect the O-term in (5.2) for n 2 N, the requirement 
(5.1) may be simplified to 


—p 
e-2V 202 <_< 10 


a ry erV2Nz , 


or even to 


(5.3) e-tvine < 19 _ 


having made the right-hand bound smaller. Inequality (5.3) yields 


ss (yes + pln 10 + es) 
= 2./ax ’ 


which gives us the desired estimate of ». 





Oak Ridge National Laboratory 
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Expansion of Hypergeometric Functions in 
Series of Other Hypergeometric Functions 


By Yudell L. Luke and Richard L. Coleman 


Abstract. In a previous paper [1] one of us developed an expansion for the con- 
fluent hypergeometric function in series of Bessel functions. A different expansion 
of the same kind given by Buchholz [2] was also studied. Since publication of [1], 
it was found that Rice [3] has also developed an expansion of this type, and yet a 
fourth expansion of this kind can be deduced from some recent work by Alavi and 
Wells [4]. In this note, we first deduce a multiplication formula for the Gaussian 
hypergeometric function which generalizes a statement of Chaundy, see (11), page 
187 of [5], and includes a multiplication theorem for the confluent hypergeometric 
functions due to Erdélyi, see (7), page 283 of [5]. Our principal result is specialized 
to give an expansion of the confluent hypergeometric function in series of Bessel 
functions which includes the four above as special cases. With the aid of the Laplace 
transform, the latter result is used to derive an expansion of the Gaussian hyper- 
geometric function in series of functions of the same kind with changed argument. 
This is advantageous since, throughout most of the unit disc, the change in argu- 
ment leads to more rapidly converging series. For special values of the parameters, 
the expansion degenerates into known quadratic transformations. 


1. A Multiplication Theorem for the Gaussian Hypergeometric Function. We 
first prove 


, (ab), \_ S(—)*(a)(8)e 2 ie ) 
ar ( Je) = > Ha+h “"\ «ape |* 


, (atk pBt+k! 
For notation and definitions, see Erdélyi [5], Chapters 2 and 4. Note that the param- 


eters a, 8 and y are free, provided that the resulting expression is meaningful. Let 
dm be the coefficient of \” on the right-hand side of (1.1). Then 


aan (= )'2'(a)e(B)a(—k)mCy + k)m(@)m(b)m ¥ +k, B+k :) 


(1.1) 








es, ky + k)u(a)m(B)m(C)m m! y+1+2k 
_ (@)m(b)m 2” $ (ae + m)n(B + m)n 2” S (—)'T(y + 2m + k) 
~ a oo n! to «OKT (y + 2m + 2k) 


(1.2) 
_ My +1 + 2m + 2k)T(n + 1) 
(iy +1+2m+n+ k(n + 1 — k) 
- (a) m(0) m 2” = (a + m) (6 T mM)» z" B 
(c)mm! ‘tao mi(y +1 + 2m), vty 
Received September 30, 1960. This research was supported jointly by the United States Air 
Force through the Air Force Office of Scientific Research of the Air Research and Development 


Command and by the United States Navy through the Applied Mathematics Laboratory of the 
David W. Taylor Model Basin. 
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where 


a —n,y+2m | 
SE) ert en ) 
(1.3) 





m _2n F + 
y+1+2m+n**\ y+ 
Now using the fact, see (46), page 104 of [5], 


F, . B 1) _ r(c)r(c — A — B) 
(14) ? Cc |\'J TC—ArC —B)’ 


C #0,-1,-2,---, R(C—A-—B)>0, 


we see that 
il f n=O 
(1.5) B, = 
0H a> O. 
Thus (1.2) becomes 
» a) m (0) m2 ¢ 
(1.6) Lgl 
' "+ 


which proves (1.1). Chaundy’s statement is (1.1) with A = 
In (1.1), set b = B, replace z by z/b and let b—> ~. Then 
= (— )*(a)e 2 
taa%) = > + 
(a, c; dz) 2 Ely + ba 
(1.7) aa 
- 3F (~ indie 2) 4), a4 +ky+1 + 2k;z). 


a,c 


Erdélyi’s expansion theorem is (1.7) with a = a. Again in (1.1), put b = c, replace 
d by A/a and let a > ~. Then 


ds (—)*(a)(Bee , (—k,k+y7 ,fat+tkB+k 
(1.8) e =>) Get by. Ps( a, B SY litle get z). 


It can be noted that expansion formulas of the type (1.1) for the G-function (a 
generalization of hypergeometric functions) have been studied in a series of papers 
by Meijer [6]. However, neither (1.1) nor (1.7)—(1.8) can be deduced from his work. 


2. Expansion of the Confluent Hypergeometric Function in Series of Bessel 
Functions. In (1.7), put y + 1 = 2a and note that 
(2.1) (a + k; 2a + 2k; 2) = Pla +k + 3) (4/2) Taye apa(2/2) 
Then 





(a, ¢; dz) = oe Is(2/2) 
(3)32 
(2.2) id , ‘ 
+7 . 2 > inn ne. oe sh R,(a, c, 6; X)In4s(2z/2) 





fo 


ace 


1 (a 
pers 
ork. 


ssel 


2/2) 
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where a — 4 = 6 and 


F . y —k,k + 26, a | 
(2.3) R,(a, Cc, 6; Xr) = Fs ( 4 + 3, ¢ a). 


Assume \ = | and write R,(a, c, 5,1) = R.(a, c, 6). If & = 0, we get the result in 
[1], and 6 = a — 3 yields the Buchholz expansion. Another proof of the latter has 
been given by Slater [7]. 6 = 3 and 6 = c — } give the expansions in [3] and [4], 
respectively. The 3/2 given in [4] is not of the form (2.3), but may be reduced to it 
using known transformation properties of ;/2’s (see Bailey [8]). 

Properties of the R,’s for \ = 1 are next of interest. These are recorded without 
proof as the argument is much akin to that developed in [1] for the case 6 = 0. 


(2.4) R.(a, c, 8) = (—)*Ry(e — a,c, 8); Ri(a, a, 8) = (—)*. 
(k a cjy(k a 26) Resi(a, Cc, 5) 


(2.5) 
= 2(c — 2a)(k + 6)R,(a, c, 6) + k(k + 26 — c)Ry_(a, c, 8).* 


. (4),.(46 — a + 4), 
(26) R,(a, 2a,6) = 0, k odd; Ry(a, 2a,5) = ~2— co . 
es . G + ald + ds 


\ (c — 2a), 
27 R,(a, c,c — a — 3) = - 
( ' , (c)x 
ki 
(2.8) R.(a, c,a — 4) = S—) as — Che 
(c), 
If 6 = 3, then 
(2.9) Erf z = [ on dt = og! STi (27/2). 
0 k=0 


The convergence of (2.9) is inferior to that of (2.7) of [1]. 


3. Expansion of the Gaussian Hypergeometric Function in Series of Functions 
of the Same Kind with Changed Argument. Now combine the known Laplace 
transforms [9] 


ne 2 b)(p — 4)" o—-s 
P17, (4) dt = Te + MP @) op (1 — 0,050 +1; 
(3.1) I . . I,{t) dt T(v + 1)¢ oF; 1 ); but 3 > . 


q=(p—1)"*, R(b)>0, R(p)>1, 


_ . r'(b) ) 
(p+1)t,b—1 " i. 
e t @(a;c; 2t) dt = —— F («, b¢ ), 
(3.2) [ (p+ 1)>? 1 D o> 


R(b) > 0, R(p) >0 





* A proof after the manner of (2.5) shows that R,(a, c, 5; \) satisfies a five-term recurrence 
formula. The expression does not seem to be of general interest and so is omitted. 
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with (2.2) for’ = 1 and put z(p + 1) = 2. Then for |z| < 1, 


w2™ & (—)*(b), (28), 1—w\ 
— Pa co 8 é ie i 
F(a, b;¢52) = es 2 Ke), Pee 8) ( + “) 





(3.3) 
x oF; (1 ~b44d-68+041,- 22) 
4w 


where w = (1 — z)”” is real and positive if z is real and 0 < z < 1. Equation (3.3) 


is convenient since the change in argument leads to more rapidly convergent series 


throughout most of the unit disc. To see this, note that |z| = re =f for all 

|z| < 1 and equality holds only for z = 0 and z = 1. Also, if z is real, | z| = 
1/2 

(1 — w)?/4w = v whenever z S i am 0.96924. If z = pe”, then by 


32 
numerical computation, we find that | z| = | v | holds for p = 0.97, 0.98, 0.99 and 
1.0 if @ = 0.0074, 0.0260, 0.0331 and 0.0365, respectively. Another advantage is that 
6 is a free parameter. Thus, if 6 = b or 6 = b — 1, the 2F; on the right is unity. If 
also the parameters a, b and c are such that R, is a product of gamma functions as 
in (2.6)—(2.8), then the right-hand side of (3.3) is a hypergeometric function and 
we can recover known quadratic transformations. For example, if 6 = 6b, and 
c = 2a, then 


(3.4) F(a, b;2a;2) = (—2_) ar, ( Fat +ale 
oO. 2f\\ a, 0, 24,2) = i+w of) —O45,55 +4; l+w 


and from the latter, one can readily deduce the Gauss transformation for the com- 
plete elliptic integral of the first kind. 

If a = 3,b = 1,c = $ and 6 = 0, (3.3) yields expansions for are tan a and 
log (1 + a) which are special cases of Chebyshev expansions previously given by 
Luke [10]. 

We now obtain a useful expansion for the complete elliptic integral of the first 
kind K(k) when the modulus is near unity. It is known that 


K(k) + 


ain 


:(@,) 
(log k’)K(k’) = 30 \3 “f 


(3.5) mm | om! 


f | . > 
“vim + 1 ) _ vy (m + 5) (k’)™, (k’)* =li1— ce, 
which is the same as 


(3.6) nie E (« ~ anaes ey) + * (log 2)K(K’). 
da , a=1 T 


To evaluate the partial derivative, use (3.3) with 6 = 0 and b = 3 where first a is 
replaced by a — 3 and then ¢ is replaced by a. It follows that 


ee . 
(3.7) 2s ( - },a,0)} =i ~ 5° k>0 
da 2 — | = 





m- 


ind 
by 


rst 











EXPANSION OF HYPERGEOMETRIC FUNCTIONS 237 


and so 


2 4 —1/2 “G 5). i—& 
(38) K(k) = (tog *) K(k’) + 2k el — ( od 


m=1 !m 1+k 


a 8 (1 — k) 
‘hh ; tere ope 


Midwest Research Institute 
Kansas City, Missouri 
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Solution of the von Mises Boundary Layer 
Equation Using a High-Speed Computer 


By A. R. Mitchell 


1. Introduction. In the problem of incompressible boundary layer flow along a 
semi-infinite flat plate against an adverse pressure gradient, the non-dimensional 
free stream velocity is taken to be 


um =-1—2Z 


where the plate lies along the positive part of the x-axis. The von Mises equation 
for incompressible flow is, in the usual non-dimensional notation, 





dz az 
(1) an Op 
where 
z2=u —w. 
On the plate, = = —2(1 — x) and uw is zero and so from (1), “4 is infinite. This 


unpleasant singularity on the plate caused Gértler [1] to abandon finite difference 
methods of solution of the von Mises equation and it looked as if the comparative 
simplicity of (1) could never be utilized from the point of view of numerical solution. 
Thomson and the present author [2], however, showed that finite difference solutions 
are possible using the von Mises equation and obtained the expansion 


‘ 2 (dz 7 (az\ 2 2 
9 ‘a 1/2 ie pra os ae ~ 3/2 ' 2 Ayia 
(2) —— ~“s (=) ¥ — ipa (=) . oe 


for the velocity in the vicinity of the plate, where a depends only on x. This expres- 
sion for u incorporates the conditions of compatibility at the plate and provides a 
means of obtaining a value of the skin friction on the plate from the computed values 
of the velocity in the boundary layer. The skin friction is given by the value of 





~ xi on the plate, which from (2) is equal to 4a’. 

The calculation in [2], however, was carried out on a desk computer and its ex- 
tent was limited, particularly in the vicinity of the plate. In the present paper it is 
hoped to remedy this defect by using a finite difference form of the von Mises equa- 
tion which is particularly suitable for exploring the region near the plate, and 
carrying out the calculation on an automatic computer. 


2. The Finite Difference Equation. The four-point explicit difference replace- 
ment of (1) with y = rfand x = sh is 


ere = A2rst,s + Bz,.s + ening 
Received November 21, 1960; revised February 10, 1961 
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where A = C = 6u,,,, and B = 1 — 26u,,,, with 6 = h/¢’, where h and ¢ are the 
non-dimensional mesh lengths in the z- and y-directions respectively. Since 
A+B+C=1,A>0,C > 0, it follows from Richtmyer [3] that z is bounded 
provided B > 0, which leads to 


ieee 


2u,,. 


This is a very satisfactory condition, since it means that in the region next to the 
plate where u is small, h can be chosen correspondingly large for a given value of ¢. 
From the point of view of computation, however, it is much better to keep h con- 
stant, and so ¢ can be reduced in the vicinity of the plate. This enables the region 
next to the plate to be examined in greater detail than the rest of the field. 

With this in view, a distribution of nodes in the y-direction is chosen according 
to the formula 


(3) vy, = $r(r+1)% (r = 0, 1, 2, 3, --- DD). 
In the x-direction, the nodes are given by 

Z, = oh (e = 0, 1, 2, 3, ---). 
The four-point forward difference replacement of (1) is now 


Terie + (r + 1)era.. — 2(r + 3 )zr.0 
rr + A) + 1) 


where 6 = h/f:' and (4) can be rearranged to give 


(4) 2r,s+1 = Zr. + ta 





2r,s+1 _ De2r41,6 + Ez,,s + Pants 





where 
ee 
(r+ 3) 4+ 1) 
; 26u,.s 
eel ore 
bu 
ee. r,s > 0 
rr + 4) 
andD+H+F=14.If 
i. rir. + 1) 
(5) 6 <- Dur. 


it follows that E > 0, and so the theoretical solution of the difference equation (4) 
is bounded. However, in view of the non-linear nature of the difference equation, 
condition (5) does not necessarily guarantee the boundedness of the numerical solu- 
tion of (4), and a watchful eye’should be kept on the calculation in case any instability 
arises. 
The principal part of the truncation error in (4) is given by 
. h oz = bly” oe . 
2 @& 3 oy’ 
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In the vicinity of the plate, from (2), 


az 1 3/2 
—=-(U-2 a 
ts ( yy 
and so the truncation error at nodes near the plate becomes approximately 
—}#(1 r x)h. 


3 
Thus, although = and higher derivatives may be very large in the neighborhood of 


the plate, the truncation error can in fact be small if h is small. 


3. The Calculation. In the calculation carried out in the present paper, the 
starting values of z at « = 0.05 are obtained from Howarth [4] after a certain 
amount of computation. The values to eight places of decimals are given in Table I. 
The mesh lengths are h = 0.0005 and f = 0.01, and so 6 = 5. It can be seen from 
the last column of Table I that the inequality (5) is satisfied for all r at the start 
of the computation. Since the values of u decrease subsequently at nodes in the 
same row, there is no danger of condition (5) being violated anywhere in the field. 

The calculation was carried out explicitly using (4) on a DaTaTRON 205 high- 
speed computer, and values of z obtained were rounded off after eight places of 
decimals. The machine was allowed to run until a value of z was obtained on 
¥v = 0.01 which exceeded the boundary value on y = 0 at the same station x. The 
machine was then confronted with the task of taking the square root of a negative 
quantity, and so ceased to print out. This actually occurred at « = 0.137. Thus, 























TABLE I 

a ard ; een 
0 0 0 0.90250000 

1 0.01 0.1480 0.88059600 | 6.76 
2 0.03 | 0.2620 0.83385600 | 11.45 
3 0.06 | 0.3726 0.76366924 16.10 
4 0.10 | 0.4846 0.66766284 | 20.64 
5 0.15 | 0.5910 0.55321900 sis 25.38 
6 0.21 | 0.6845 0.43395975 | 30.68 
7 0.28 | 0.7673 0.31375071 | 36.49 
8 0.36 | 0.8337 0.20744431 | 43.18 
9 0.45 | 0.8831 0.12263439 50.96 
10 0.55 | 0.9165 0.06252775 | 60.01 
11 0.66 | 0.9355 0.02733975 | 70.55 
12 0.78 0.9452 0.00909696 | 82.52 
13 0.91 0.9487 0.00246831 | 95.92 
14 1.05 | 0.9499 0.00018999 | 110.54 
15 1.20 | 0.9500 0.00000000 | 126.32 
16 1.36 | 0.9500 | 0.00000000 143.16 
17 1.53 0.9500 |  0.00000000 161.05 
18 1.71 ' 0.9500 | 0.00000000 180.00 
19 1.90 0.9500 | 0.00000000 200.00 
20 2.10 | 0.9500 |  0.00000000 221.05 
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TaBLeE II 
=z Zz Mm 
0.0500 | 0.88059600 0.1480 
0.0600 | 0.86538596 0.1349 
0.0700 0.84951282 0.1244 
0.0800 0.83344382 0.1125 
0.0900 0.81730237 0.1039 
0.1000 0.80118136 0.0939 
0.1100 0.78516342 0.0833 
0.1200 0.76935383 0.0713 
0.1250 0.76158997 0.0635 
0.1300 0.75401304 0.0537 
0.1350 0.74695089 0.0357 
0.1355 0.74632746 0.0321 
0.1360 0.74574286 0.0274 
0.1365 0.74521891 0.0203 





values of z were obtained in the field from x = 0.05 to x = 0.1365 at intervals of 
0.0005 for twenty values of y between y = 0 and y = 2.10. The values of z, to- 
gether with the corresponding values of u obtained at y = 0.01, are shown in Table 
II for several stations x. The actual running time of the machine was determined 
by the time required to print out the results. When only the values of z at y = 0 
and y = 0.01 were printed out, the machine took only four minutes for the entire 
run. This gives a good indication of the simplicity of the present scheme for solving 
the boundary layer equation. 


4. Comparison of Results. Accuracy must not be sacrificed for speed and 
simplicity of calculation, however, so the results are now compared with those ob- 
tained by Leigh [5], who solved the boundary layer equations numerically with u as 
a function of x and y, using the approximation of Hartree and Womersley [6]. It is 
sufficient to say that the results in the present calculation at x = 0.12 agree with 
those of Leigh at x = 0.1198 to within 0.5 per cent. For example, at y = 0.01 
x = 0.12, the node in the present calculation where the disagreement is likely to be 
greatest, the value of u is 0.0713, compared with a value of u between 0.071 and 
0.072 from Leigh at Y = 0.01, = 0.1198. It is unwise to quote Leigh’s results more 
accurately at the nodes of the present calculation, as numerical integration and 
interpolation are necessary to obtain them from his original calculation. The close 
agreement is very encouraging considering the vastly different natures of the two 
calculations. Leigh started his calculation at « = 0.10 with values taken from 
Hartree [7], whereas the present calculation commenced at x = 0.05 with values 
taken from Howarth. Also, Leigh solved a set of simultaneous linear equations at 


each step x, whereas the present calculation involves a simple explicit calculation at 
each node. 


5. The Singular Solution. Returning to the calculation in this paper, it is inter- 
esting to examine the solution in the neighborhood of the breakdown station 2» which 


is somewhere between x = 0.1365 and x = 0.1370. By plotting log| — (*) | 
4 ¥=0 
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against log (x, — x) for different values of x, , it is found that 


[- 3 (G5),.]<» — = 
3 av) yno| **” M 


where the corresponding values of 2, and q are 


Xp 0.1369 0.1368 0.1367 0.1366 
q 0.73 0.67 0.60 0.52 


The values of x considered for each x, are 0.1365, 0.1360, 0.1355, 0.1350, and at 
each of these stations (*2) is given by 


y=0 


(*) _ 0.1 — 2y—0 
Ay / y= 0.01 : 


Despite the fact that the above result is based on a breakdown point in a finite 
difference calculation, there is a distinct resemblance between it and the result 


\ 
(2 a(x, — x)! 
dy] y=0 


obtained by Goldstein [8] from the asymptotic solution of the differential equation 
valid in the neighborhood of the separation point, where x, is the separation point, 


and (2) is the skin friction in the physical (x, y) plane. Consequently, it is felt 
J/ y=0 


that further finite difference calculations using the von Mises variables may go a 
long way towards determining the nature of the solution of the boundary layer 
equation in the neighborhood of the separation point. 
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The Numerical Solution of Coupled Integro- 
Differential Equations 


By M. M. Pennell and L. M. Delves 


1. Introduction. The purpose of this paper is to report on a method for the nu- 
merical solution of simultaneous integro-differential equations of the form 


[ X Ge, 99) dy = Palade) 
0 


n=0 n=0 


[ & Gea, 072) do = Y PaCr)g(r) 
0 n=0 n=0 

where the asymptotic forms of f(x) and g(r) are known; and where ,k, and 2k, are 
the kernels. The method, which is reasonably simple to program for a computer, has 
been tested on a problem arising in nuclear physics and yielded reasonable results. 


2. Motivation. A great deal of effort has gone into evolving satisfactory numeri- 
cal methods of solutions of differential equations. Much less has been done in the 
field of integral equations; Fox and Goodwin [1] have given methods for solutions 
of equations of the form 


b 
(2.1) / k(x, y)f(y) dy = p(x) + f(x) 


where p(x) and k(z, y) are known functions. Their method is to write the left-hand 
side of equation (2.1) in its finite difference form 


m 


b=a+mh 
(2.2) [Ka wy) dy = XS Ay k(x, nh)flnh) + 


n=0 

where A, are the coefficients and , the truncation error of the particular quadrature 
chosen. The substitution of equation (2.2) into equation (2.1) for z = a,a + h, 
a+ 2h, ---,a-+ mh, gives m + 1 simultaneous equations in x, and the m + 1 
unknowns f(a), f(a + h), --- , f(a + mh). The truncation errors x, may be made 
negligible by suitable choice of step-length or treated by the iterative scheme given 
by Fox and Goodwin; the equations may then be solved by any of the standard 
methods. 

In the same paper they note that this method extends itself immediately to 
integro-differential equations of the following form 


b m 
(23) [ Hz, WI) dy = P(x) + fle) +O pal2f C2) 


where f‘”’(x) is the nth derivative of f(x), and k(x, y) and p,(«) are known fune- 
tions. The only additional substitution needed is the appropriate finite difference 
approximation for the derivatives. The treatment of the end points depends upon 
the boundary conditions, but in general causes no trouble. 


Received October 18, 1960. This work was supported by the U. 8. Atomic Energy Commis- 
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There arise in nuclear physics equations of the following form: 


(2.4) [ Dd (kala, g)g™(r)) dy = ay iP,(x)f\ (x) 
0 n=) n= 
(2.5) [ &X Gebel, oF (@)) do = YE aPalr)9 (2) 
0 n= n=! 
where x, y and 7, @ are related coordinates 
x = a(r, 6) 
(2.6) 
y = p(r, 6). 


The problem is to compute a table of values of f(x) and g(r) for x = r = a, 
r=r=at+h,--:,x =r=a-+ mh (these values of x and r will be referred to 
as the mesh points). 

Equations (2.4) and (2.5) differ from equation (2.3) in the following significant 
respects: 

a) The appearance of the derivatives of the unknown function in the integrand 

b) The infinite interval of integration; that is, the equations are singular in the 
sense of Fox and Goodwin 

c) The appearance of two naturally occurring coordinate systems; a finite 
difference mesh at even intervals in one coordinate system does not in general repre- 
sent even intervals in the other. This reflects the fact that (2.4) and (2.5) are two- 
dimensional equations, although only of a restricted type. 


3. Method. The method of Fox and Goodwin can nonetheless be used, but with 
appropriate modifications, to take care of the above difficulties: 

a) The derivative under the integral sign may be replaced by an appropriate 
finite difference approximation; for example, the second difference might be used 
for the second derivative. 

b) The infinite integrals can be made finite if there is a point beyond which 
the contributions of the integrand are negligible. However, in general, this will in- 
volve values of x and r outside the range in which the solution is desired. 

c) Any two of the four coordinates determine the other two by relation (2.6). 
For example, in equation (2.4) r takes on the values a, a + h, --- , a + mh; but 
for a given r it is necessary to integrate over y. Each combination of r and y gives 
by equation (2.6) a value of x which may or may not be a mesh point. Two cases 
naturally follow: 


a) a<x<a+mh but +#a+ith imsidethemesh 7 =0,1,2,---,m 
b) xz<a or a+mh<<x_ outside the mesh. 


The first difficulty may be taken care of by interpolation. If, for example, linear 
interpolation is chosen and z lies between the mesh points x; and 24: 


f(x) = f(xi)q + f(xiss)P 
ie - Be 


=l1-—- = —— <2 < 
Pp q Pp ine +1 
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Similarly, after substitution of the second difference, 
f(x) = pf(xis2) + (3q — 2)f (igs) + (3p — 2)f(a) + af(tin) 44 <2 < 24. 


Alternatively, if x lies beyond the end values of the mesh, it is possible to em- 
ploy the (assumed known) asymptotic forms of f(z) and g(r). Here the assumption 
is made that the last two mesh points tr = r = a + (m — 1)handz = r = a+ mh 


are in the range of validity of the asymptotic forms. Suppose, for example, the 
asymptotic form of f(z) is known to be 


(3.1) f(x) — Agi(z) + Bdn(zx) r— © 


where ¢; and ¢ are known functions; and A and B are unknown constants to be de- 
termined a posteriori from the solution of the problem. Thus by our assumption 


f(a + (m — 1)h) = Agi(a + (m — 1)h) + Bola + (m — 1)h) 
f(a + mh) = Adi(a + mh) + Boo(a + mh) 


A and B can thus be solved for in terms of the two unknown end values 


f(a + (m — 1)h), f(a + mh); and these expressions can in turn be substituted 
back into equation (3.1) giving 





f(z)— ; 
pes — (m — 1)h)¢.(a + (m — 1)h) 
¢:(a + mh) ¢2(a + mh) 


f(a+ (m — 1)h)o.(a+ (m—1)h) | 
f(a+mh) d(atmr) |%2) + 


o:(a +(m — 1)h)f(a+(m — 1)h) 
¢i(a+mh) fla+mh) | (2)) 


( 








1 
= _ , h 1 
i +o nO ea - te eet 


¢:(a + mh) ¢2(a + mh) 


1 
~ 0+ RE + Ke a aie oe i) 
é:(a + mh) $:(a + mh) | 


-(oila + (m — 1)h)o2(x) — do(a + (m — 1)h)¢Gi(z)). 








Thus if x is not a mesh point, it is possible by interpolation and the use of the 
asymptotic form to express f(z) in terms of the value at mesh points. In the same 
fashion the difficulties are overcome in equation (2.5) when r takes on the values 
r=a,r=a+th,---,r=a+mh. 

The appropriate mesh size, interval of integration and cut-off point for the in- 
finite integral depend upon the problem. 


4. Example. In the discussion of neutron-deuteron scattering by the resonating 
group method including inelastic channels [3], the following equations arise: 
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SE) 5 2m e — wage) + MOI) (Gye, y) ay 
wit 15 ~9/2 2m 
(4.1) +2 if (- ee a(r)) o(y)ydy + me 


[ G,(x, y)g(r) dy = 
0 


16 1/2 3d t( x) d’g(r) 
ae r ‘t ‘Ge + G.(r, ann ¢(y) sin a cos a da + — a 





2 


sa Fae) 4. tr) + sr 9(r) [ Gi(r, «) sin? @ cos* a da = 0 
0 


x= Vir cosa y = V2rsina 


? _ oly) oR “es ae 
2B ee fom (-3) 


G, = - y [z + Vy “ts - st 














T ro xy 


eS. 
ee 
a 
| a 





fo)“ 


The physical interpretation of the various quantities is as follows: @(y) is the 
bound-state deuteron function satisfying the equation 


: . F ’ 
E aa si (1 9 exp (-4) + Ba) | o(y) = 


For large y it falls off exponentially: 
o(y) — A exp (—kay) Wha = —mE, 
The function ¢(y) and the eigenvalue Vo were calculated numerically by a 


separate program. 
Ez is the (negative) deuteron energy; m the mass of a nucleon. The potential 


, : » r 
between two nucleons is the Gaussian form — Vo exp (-5 


2 


The unknowns f(z) and g(r) represent respectively the elastic and inelastic 
scattering wave functions. 
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Two Independent Solutions (f,,9, ond f2,92) of Equations [4.1] [4.2] 
with Boundory Conditions f(O) =g(O) =O; f,(5) = g)(5) «-1 

fo(5) = -1, go(5) = +1 
fois the Solution of [4.1] Obtained by Assuming g(r) =O 


| Fie. 1. 
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x,r, fermis 
he fa eR 
i + Sitiemaescnanate 
-2- 
Two Independent Solutions (f,,g; ond f2,92)0f Equations [4.1],[4.2] 
with Boundary Conditions f(O) = g(O) =O; f, (5) = g, (5) =-1 
. = fo(5) =-1, g2(5) = +1 
fo is the Solution of [4.1] Obtained by Assuming g(r) #0 
a) Fie. 2. 
‘ial The equations are of the form of equations (2.4) and (2.5), and have been 
solved by the methods described here. The asymptotic forms of f(z) and g(r) are 
stic f(x) = Asin (k,x) + B cos (kar) Qk, = 3m(E — Ea) rt —- & 
=) 


(4.4) 


g(r) = Cr? J2(kr) + Dr? N2(kr) Wk? =2mE r— 
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where J,,(z) and N,,(z) are the Bessel functions of the first and second kind of order 
n as defined in Jahnke and Emde [2]. 

The application of the above techniques for solving the equations is straight- 
forward. The integrals of infinite range converge rather slowly, so that the back 
interpolation (equation (3.2)) has to be applied over a wide range of r; however, 
this can be avoided in this example by noting that for r such that equation (4.4) 
is valid 


2 
(ns jt ty 2” Gs) g(r) = 0 


dr? 49/2 h? 


since this reduces to the Bessel equation satisfied by g(r). The second and third 
integrals in equation (4.1) can then be cut off at this point, which is long before the 
individual terms are negligible. Moreover, the first integral converges very rapidly 
and is negligible even before this cutoff point is reached. The equations have been 
solved in particular cases both with and without taking note of this point; the results 
were comparable but the computing time was reduced from fifteen to three minutes. 
Difference corrections were neglected; the mesh size was varied until stability was 
reached with the accuracy required (~3% ). Figures 1 and 2 give a plot of the results 
of runs for E = 0, E = 7; the constants used were (lengths in units of 10-” em; 
energies in Mev) 


To = 1332 V> = 86.674 Ea = —2.226 = = .0481933 


h(y) = 2 h(x) = 25=h(r) E=5. 


The asymptotic form was used for r, « = 4.75. Other runs were made for E = 
0.5, 1, 3, 5; for E = 0, the asymptotic form for g(r) is modified. Also plotted in 
Figures 1 and 2 are the results obtained by solving equation (4.1) for f(x) setting 
g(r) = 0 (physically, neglecting the effect of the inelastic channel upon the elastic 
channel). For further detail of the numerical results see [3]. Higher values of E were 
not used because the end point of the mesh used (r = 5) is close to a zero of g(r) 
for E of the order of 10. 

The problem was coded in SAP for the IBM 704 at M.I.T. Coding was straight- 
forward; the most laborious task was debugging—a single suitable hand check took 
several days to complete. The program was approximately 2,000 instructions long. 
Running time for the coefficients from (4.1) was from 2.5 to 7 minutes depending 
upon mesh size; 2 minutes for the coefficients from equation (4.2), and one minute 
for the matrix solution. 


5. Acknowledgements. The authors wish to acknowledge the help and advice of 
Miss E. J. Campbell of the Joint Computing Group, M.I.T. Mr. Delves is indebted 
to the Shell Petroleum Company for the award of a Shell Scholarship and to Pro- 
fessor V. Weisskopf for making possible his stay at M.I.T. This work was done in 
part at the M.I.T. Computation Center and in part at the M.I.T. Laboratory for 





ce of 
bted 
Pro- 
ne in 
y for 





NUMERICAL SOLUTION OF COUPLED INTEGRO-DIFFERENTIAL EQUATIONS 249 


Nuclear Science, Massachusetts Institute of Technology, Cambridge, Massa- 
chusetts and Clarendon Laboratory, Oxford, England. 


Clarendon Laboratory 
Oxford, England, and 


Joint Computing Group 
Massachusetts Institute of Technology 
Cambridge, Massachusetts 


1. L. Fox & E. T. Goopwin, “‘The numerical solution of non-singular linear integral 
equations,’ Philos. Trans. Roy. Soc. London, Series A, v. 245, 1953, p. 501-534. 

2. JAHNKE & Empe, Funktionentafeln; Teubner, Berlin, 1933. 

3. L. M. Detves, “Elastic and inelastic neutron-deuteron scattering near the inelastic 
ae ye submitted to Nuclear Physics. 

4. L.F 


. Fox, The Numerical Solution of Two-point Boundary Problems in Ordinary Differential 
Equations, Oxford, 1957. 








A Computational Approach to the Four-Color 
Problem 


By Hidehiko Yamabe and David Pope 


Introduction. The problem of coloring a map on a sphere with four colors was 
shown by Heawood [3] in 1898 to be equivalent to the solution of a certain set of 
congruences, which we will call the assignment problem in this paper. The assign- 
ment problem is shown to be amenable to a simple computational searching routine, 
which can be easily coded for a binary automatic digital computer. A brief descrip- 
tion of such a code for the UNIVAC Scientific 1103 computer is given, with two 
examples. 

This work was completed just before Professor Yamabe’s untimely death on 
November 20, 1960. It was hoped that the experimental results found by the com- 
puting machine might lead to an induction hypothesis, in particular, one involving 
the complete set of solutions to the assignment problem. 


1. The Assignment Problem. By a “map” we will mean here a polygonal parti- 
tion of the 2-sphere into a finite number of open simply-connected domains (‘‘coun- 
tries’). Each country will be bounded by a polygon with a finite number of edges 
and vertices. Two countries will be called adjacent if they have at least one edge in 
common. It is easily shown to be sufficient to look at the case when any two adjacent 
countries have only one edge in common. By a vertex we mean a point where at 
least three edges meet. By a coloring of the map we mean an assignment of one 
color to each country such that no two adjacent countries have the same color. The 
coloring problem is easily reduced to the case of a normal map, that is, where each 
vertex meets exactly three edges (Cf. [4], p. 298). In the rest of this paper we shall 
consider only normal maps. 


We now define an assignment of the map as a function ¢ of the vertices of the 
map, giving each vertex v of the map the value ¢(v) = +1 or —1. Then we define 
a proper assignment of the map as an assignment such that for each country, the sum 
of the values ¢(v) taken over the vertices on the boundary of that country is congru- 
ent to zero modulo 3. We then have the following theorem, due to Heawood [3]. 

THEOREM 1. A normal map can be colored with four colors if and only if there exists 
a proper assignment of the map. 

To derive the four-coloring of a map from a solution to the assignment problem, 
we first label the edges of the map a, b, or c such that at each vertex v three differ- 
ent letters meet, and the ordering of the edges defined by a — b — ¢ is a rotation 
in the positive direction around the vertex v if ¢(v) = +1, and a rotation in the 
negative direction if ¢(v) = —1. This labeling is consistent if and only if the assign- 
ment is a proper one. 

Given such a labeling of the edges, the closed curves formed by the edges labeled 
alternately a and b separate the sphere into two disjoint (but not necessarily simply 
connected) regions, to which we may assign the numbers 0 and 1. The closed curves 
formed by the edges a and c also separate the sphere into two regions, to which we 
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assign the numbers 0 and 2. The sum of these two numbers gives us the four-color- 
ing of the map expressed as the “colors” 0, 1, 2, 3. (For a further discussion and 
proof, see [1] p. 226 and [2] p. 13.) 


2. Description of the Code. The solutions to the assignment problem for the 
vertices of a given map are found by a searching routine, which was coded on the 
unIvac Scientific 1103 computer at the University of Minnesota Scientific Computa- 
tion Laboratory. It is of interest to note that for this problem a binary computer 
is more efficient than a binary-coded decimal computer. 

Each vertex of the map is represented in the computer by one of the 36 bits of 
an 1103 word; thus the present code is restricted to 36 or fewer vertices. Each coun- 
try is then defined in one word, containing a binary “1” in each of the positions 
representing a vertex lying on its boundary, all the other binary digits being zero. 
The vertices are numbered from 1 to 36, and vertex j is represented by the jth binary 
digit from the left. 

An assignment of +1 or —1 to each vertex of the map is also given in a single 
word, called A in the code, with +1 represented by a binary one, and —1 repre- 
sented by a binary zero in the bit corresponding to the given vertex. Thus every 
possible assignment can be formed in A by successively adding one to the least sig- 
nificant digit, running A from zero to 2° — 1. (The complement form used for nega- 
tive numbers is circumvented in this code.) 

In practice, of course, this method of assignment generation would be highly in- 
efficient. However, we may eliminate many assignments at once by a proper ar- 
rangement of the countries and vertices. To describe this arrangement, let v(k) be 
the vertex with highest index lying on the boundary of country k. We then rear- 
range the ordering of the countries so that if kj < k. , then v(k,) S v(ke). We then 
check the given assignment, beginning with country 1. If the assignment does not 
satisfy the congruence for country 1, we may add 2°" to A instead of 1, since the 


4 


Fic. 1.—Map One 
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Fic. 2.—A Map With 20 Countries 


bits bys.) --- bo of A do not affect this congruence. If this assignment does satisfy 
the congruence for country 1, we proceed to country 2, and so on. In this way many 
possibilities are eliminated at once for the countries with low index. The assign- 
ments satisfying the congruence for all countries are then printed out. Although the 
code will print out all the solutions to the assignment problem, clearly only half of 
these need to be computed as the other half are the reflection of these solutions ob- 
tained by interchanging + and —. 


3. Examples. The simplest example computed was Map One, given in Figure 1. 
The total computing time for this, including output, was Jess than two seconds. The 
input data describing this map are given below, in octal. 


Map One Inpur 


Address Contents Explanation 


00201 74 00000 00000 

00202 63 00000 00000 

00203 52 40000 00000 Code words describing countries 
00204 03 60000 00000 1 to 6 

00205 25 20000 00000 


00206 


14 60000 00000 





The output of the code gives the eight solutions to the assignment problem, as 
follows: 


Map One Ovurput 
— SN PP PS SK = 
- -+--—- + Ft -—- + - 
—-— + tre -- t+ + 
—-_ + +t =F - -— + 
+ - —- + -—- + + - 
+ - -—- + + - - + 
+ - + - - + - + 
+ + - - - —- + + 
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Note that the last four solutions are the images of the first four under the basic re- 
flection, in addition to the other symmetries generated from the symmetries of the 
map. 

The second example, given in Figure 2, has 20 countries and 36 vertices. The 
total number of solutions to the assignment problem is 146, or 73 if we discard the 
half given by the basic reflection. The first solution discovered was —-—— ——— 
—t+t+ +4 33 eee Ft Kt OH St - SH 
———, or 003 500 524 440 in the octal abbreviation. The total computation time 
for the 73 solutions was about 25 minutes. 
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Some Applications of High-Speed Computers to 
the Case n = 2 of Algebraic Cryptography 


By Jack Levine 


1. Introduction. In 1929 Hill [1] proposed the use of simultaneous linear con- 
gruences as a method of encipherment (see also [2], [3]). If the number, n, of such 
congruences be 5 or more this results in a cryptographic system of unusual security. 
In this article it is shown that high-speed computers can be used in the problem of 
the decipherment of the simplest case n = 2. We give first a brief description of the 
system using this value of n. 

The 26 letters of the alphabet are assigned numerical values according to some 
arrangement of the numbers 0, 1, 2, --- , 24, 25. For example: 


a 2G 2 4 ¢ 2 #@ 4.8 3B ££ 2 OO P CG 

72a 8 47 6 8 Pan H HSH 2 Mw 8B 
sos: C VY FP 2 YF 

1 2 203 2 5 8 2 10 


(1.1) 


A 2 x 2 involutory matrix, mod 26, is selected to form the congruences 


C, = aP, + bP,, 
C2 = cP, + dP, 


a b :.. ms 1 0 i m 
* i. Taf = & 4) mod 26. 


As an illustration we use 


4 7 
(1.4) H = E a | 


A given plain-text, say CRYPTOGRAPHY, is then divided into two-letter 
groups, CR YP TO GR AP HY; each pair of letters is selected as the PP: of (1.2), 
and (C,, C2 calculated, using the numerical equivalents of (1.1). Thus, as CR = 
P, P. = 21 1, we find, using (1.4), 


12) mod 26 


where the matrix is 


(1.3) H 


C; 


4(21) + 7(1) = 13 = N, 
Cs 


9(21) + 22(1) = 3 = U, 


so CR is enciphered by NU. The converse is also true, since matrix H is involutory. 
The complete encipherment becomes 


CR YP TO GR AP HY 
NU VN XB WJ GU LL 


(1.5) 


The decipherment, knowing the matrix, is performed in an identical manner. 
Received November 30, 1960. 
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With this in mind we may state the problem to be solved in the following way. 
Given a cipher-text, obtained as in (1.5), determine the corresponding plain-text, 
assuming as known the numerical alphabetic values as in (1.1). This amounts to 
determining the matrix H which is to be considered as the unknown quantity. As 
mentioned above, this soon becomes a complex problem with increasing n. For the 
present case of n = 2, however, it can be readily solved by machine, an IBM 650 
actually being used for this purpose. Two methods will be explained. The first 
method has the advantage that it can be applied to extremely short messages, the 
second that it can be extended with some changes to higher values of n. 


2. First Method. Since the only unknowns are the four elements a, b, c, d of the 
matrix H, the most direct method is to test in succession all such possible matrices. 
If n > 2 this becomes impractical, but if n = 2 there are only 740 matrices, a rela- 
tively small number. It is found from (1.3) that the elements must satisfy the 
conditions 
(21) @+b=1, ¢d+be=1, blat+d) =0, clat+d)=0, mod2%6 


and it is easily shown that these imply two types of matrices, 


ok, _ ja b 
rype 1. a =|° aa 


a b 
— pa 
Type 2. H [2 ’] 


with 

(2.2) a+be=1 mod 26 

in both cases. 

Type 2 contains only 8 matrices. 

Type 2 

a bee a bie 
1 0 0O 14 13 13 
1 O 18 25 0 O 
1 3 oO 25 O 18 
12 13 13 25 13 O 


The remaining 732 of Type 1 can be obtained from the basic set listed below. 


If we place 
t 
(a,b,c) = [* ae 
c —a 


then associated with (a, b, c) we have the set of eight 


(a, b, c) (—a, b, c) 
(a, c, b) (—a, c, b) 
(a, —b, —c) (—a, —b, —c) 


(a, —c, —b) (—a, —c, —b) 
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A complete list of Type 1 matrices can be exhibited by writing one of each asso- 
ciated set of 8 (or 4 or 2 in case of duplicates). Such a basic list is given below. 


Type 1. Basie List. 


a=0 a=1 a=2 a=3 a2w4 a=z5 a=6 
b c b ¢ b c b Cc b c b Cc bo ¢ 
B ud oc Lt 1 18 kL A L 2 i. a 
3 9 a. 5 15 2 9 So 2 2 14 3 2 
5 21 4 13 7 i 3 6 7 9 3 18 5 19 
2s as 6 13 > = SE 4 é | a «OE 

8 13 4 24 4 20 

10 13 5 14 5 16 

12 13 6 16 6 9 

("ec = 0, 1, 7 10 8 10 

++, 18) 8 12 11 12 
a=7 a=8 a=9 a= 10 a= 11 a= 12 a= 13 
bo ¢ b C¢ SA 2 b Cc boc b c b 6 C¢ 
1 4 1 15 1 24 1 5 1 10 | & 1 14 
2 2 3 5 2 i r 2 5 o & 2 7 
2 18 ‘oe | 3 8 9 15 2 18 5 is 2 20 
3 10 4 6 . EB 7 & a Ze 
4 14 4 19 4 9 9 13 4 10 
5 6 5 10 4 22 li # 4 23 
6 18 6 17 6 6 13 13 5 8 
7 8 8 16 6 19 6 ll 
9 12 11 14 $ 11 8 18 
10 16 10 14 9 16 
1i2 12 


Cases a = 0 and a = 13 could also be considered as Type 2, but it is convenient 
to place them here. A basic Type 2 list is given below. 


Type 2. Basic List. 


a b Cc 
1 0 0 
1 0 13 
12 13 13 


3. Machine Procedure. A card is punched for each of the 109 entries (matrices) 
in the two basic lists, and each of these in turn generates the remaining seven asso- 
ciated with it, duplicates being retained. Each matrix of a set of 8 then performs a 
deciphering operation indicated by the matrix congruence P = HC, or 


2 P, Ps Ps Pr Ps a bWC, Cs Cs Cr Cy a 
(3.1) | ] [2 ae Ce Ce Ce “A mod 26 


an een Hh Fe 
where the columns of the C-matrix represent the first five pairs of cipher letters, 
and the columns of the P-matrix the corresponding pairs of the deciphered ‘‘plain- 
text,”’ these last 10 letters being printed. When all 8 matrices of a set have been 
used, the next matrix of the basic lists is read and the entire process repeated. There 
results 8 X 109 = 872 decipherments, each of 10 letters. An inspection of these 
easily locates the correct decipherment and the corresponding matrix which is also 


Il 
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0- printed along with the decipherment. The entire procedure requires approximately 
30 minutes. 


4. Second Method. We assume here a Type 1 matrix has been used for H, all 








; Type 2 having been tested by hand. 
6 Let a cipher-text be represented as 
: (4.1) A,B AoB: «++ A;B; Ain Biss --- AmB - 
9 Suppose the plain-text contains the 3 consecutive letters P QR; and it is desired to 
1 locate their position in (4.1). The 3 letters must be divided as 
(a) PQ;R;- or 
(b) -P;Q;R;. 
Assume for case (a) the cipher-text location being tested is given by 
13 
| (4.1a) AB; AipBiss 
c 
14 PQ; R; ° 
1 Then from (1.2), in the case of a Type 1 matrix, we must have 
20 
29 
10 (4.2a) P; = aA; + bB;, (4.2¢c) Q; = cA; — aB;, 
(4.2b) A; = aP; + bQ; (42d) B;=cP;—aQ,; 
~ (4.2e) R; = aAjgi + OBin. 
1 
16 If a, b be eliminated between (4.2a, b, e) we obtain 
12 
‘ A; B; OP; 
a (4.3) Ei; =|Am Bir R;\)=0 mod 26 
| Pi QQ Ai 
In case (b), if the location be at 
A iB; A intB; +1 
(4.1b) 
*P; Q;R; 
ces) we have similarly, 
SSO- (4.4a) P; = cA; — aB;, (4.4c) Q; = dA + OB, 
ns a 
(4.4b) Biss == cQ; a aR; ; (4.4d) Bias = cQ; 7 aR; , 
(4.4e) R; — CA jas —_ aBis ° 
d 26 
Eliminating a, c from (4.4a, b, e) results in 
— iL: -* iz 
el (4.5) Fi; =| Aus Bir R; |= 0 mod 26. 
’ Q; R; Biss 
here 
hese We may state these results in the following theorem: 
also THEOREM. A necessary condition that three consecutive letters PQ;R; occur as 
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plain-text in positions (4.la) or (4.1b) of a cipher text (4.1) is that E;; = 0 or 
F;; = 0, mod 26, respectively. 

Equations (4.2a, b, c, d) can be solved to give 


_ |B Q 
= i. | Q; B; 


: _|A: Q| _|P; Ai| | 
(4.6) Da =| B,| Db = A, P,| 


’ 








In this case (4.3) is equivalent to 


(4.7) D(aA i435 + OBi4, — Rj) = 0, mod 26, 


and in addition, 
(4.8) D’(a + be — 1) = 0, mod 26. 


It follows that for this case, (a), if D is prime to 26, then a, b, c are determined 
uniquely by (4.6), with a’ + be = 1, and all equations (4.2) are satisfied for these 
values of a, b, c. 

However, if D is even, mod 26, (but not 0), several solutions of (4.6) are possible, 
and it is best to solve them mod 13 and mod 2. If a, b, c is a solution mod 13, then 
a+ 13k, b + 13k.,¢ + 13k3(k; = 0, 1) may be solutions mod 26. To pick the 
correct choices, we must solve (4.6) mod 2. There are just four solutions mod 2 
satisfying a’ + be = 1. These are 


(4.9) M, = (011), M, = (100), M; = (101), M, = (1 10), 
which when taken in conjunction with the mod 13 solutions of (4.6) will give all 
mod 26 solutions. In this case, D even, it is necessary to check (4.2e) for each such 
solution. 

If the seven numerical values under consideration in (4.la) be reduced mod 2, 


TABLE 1. (Used when E;; = 0, and D is even, mod 26) 


PQR 
ABAB 000 001 010 O11 100 101 110 111 


0000 


1,2 x Fs x x x x x 
3, 4 
0001 2,3 1,4 # x x 1 r r 
0010 1 2,3 z x x t E - a 
4 
0011 4 Es z v 1 1 x v 
3 
0100 | 2 # 2,3 x l x + rs 
0101 x 2 2, x x 1 x 4 
0110 x = = 2,3 1 2 F 2 4 
0111 x - x , x 1 4 x 
1000 Fo # 1 z 2,4 x 3 x 
1001 x x x 1 : 4 3 x 
1010 F x 1 x x 2,4 x 3 
1011 x x £ 1 4 2 F 3 
1100 F ro 4 x 3 x 2 x 
1101 F z z 4 3 s 2 1 
1110 # x x 4 S 3 2 
1111 * x 4 x x 3 v 1,2 
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then Table 1 gives the possible M; solutions of (4.2), and also all contradictions 
(these resulting from a wrong location of the PQ;R;). The entries 1, 2, 3, 4 repre- 
sent M,, M., M;, M, of (4.9) respectively, and x represents an impossible setting. 
From Table 1 we see, e.g., that the (4.1a) setting 


0 


~~ 


18 5 2: 


20 9 


_ 


converted mod 2 is 0101 010, giving M; and M; as (a bc) solutions, mod 2, whereas 
the setting 


ig 5 22 5 
20 10 4 


or 0101 001 is impossible mod 2, and therefore mod 26 also. 
Return now to the case of a (4.1b) setting. From (4.4abed) we find 


Da= | Ain: R; I, Db= Q; Ain ; 
(4.10) | Q; Ain | | Aiss Q; | 
lBus RB, | Q %#R 
D’ = i+1 j : D’ = 3 3 
, R; Biss Ais Biz, 


If D’ is prime to 26, a, b, ¢ are uniquely determined, (4.4e) is satisfied as is a’ + 
be = 1. If D’ is even mod 26 (not 0), we solve (4.10) mod 13 and mod 2. In this 
case, Table 2 gives the various mod 2 solutions for (a bc). For each (a b c) mod 26 
so obtained, (4.4e) must be tested. 


5. Machine Procedure. If an assumed trigraph P,Q;R; is actually present in the 
plain-text of (4.1) and is tested in its correct position, say A,B; Ais,Bi4, , (called 


TABLE 2. (Used when F;; = 0 and D’ is even, mod 26) 


ABAB | 000 001 010 011 100 101 110 111 





x= bo 
~ 
& 
8 | 
me 
~ 
& 
~ 


2,3 1 


bo 
_ 
oO > 


0010 | 
0011 | 
0100 | 
0101 

0110 

0111 

1000 | 2 
1001 | 
1010 
1011 
1100 | 
1101 | 


RBRReKRBRY ~ 
to 
RRRR 
w 

to 
RRRRSB 

be 


P= 
RRRReKR WY 


mMRwWHReR wr 
RRPREWRANWRHRRKRA 


1110 
1111 


RBRRWwWHRBBRRY 
RRWRPRNE RK RS Peer 
i) 
WRRARRY 
rs 
RWRRNWARKEPREKRRSRY 
Pe DR RK WR PRY 
. 
RY 
_— 
bo 
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a causal setting), one of the (a bc) solutions obtained as outlined above (assuming 
D or D’ # 0 or 13) will produce the correct decipherment. If then, a list of the 
highest frequency (English) trigraphs be tested in every cipher-text position the 
great majority of the assumed settings will be rejected. Most of those that survive 
(produce an (abc) solution) will be accidental settings, and the few remaining will 
be the causal settings. If the first few cipher pairs be deciphered using each (a b c) 
solution thus derived, the correct (a b c) matrix will be immediately evident by in- 
spection, since this will be the only one yielding English text. This correct matrix 
will be produced from the causal settings only. 
Suppose the trigraphs 


P,QiR,, --+ , PQiR;, ++: , PQeRe 


are to be tested. Such a list would include, for example, THE, AND, ING, ENT, 
HER, ION, NTH, OFT. For a given (7, 7), Eij(Fi;), (@ = 1, ---,m—1;7 = 1, 

- ,k), is evaluated and if 0 mod 26 the corresponding AB AB, PQR mod 2 entries 
in Table 1 (2) are located. If these are present (the M; solutions) then D(D’) is 
computed, and all (a bc) solutions obtained using (4.6) or (4.10) respectively (and 
using the Table 1 (2) entries in case of multiple solutions). In case D(D’) is even 
(+0) mod 26, (4.2e) or (4.4e), respectively is tested. Finally, the first five cipher 
groups are deciphered with each (a b c) solution. These decipherments, together 
with the corresponding (a 6c), are printed. If no decipherment gives plausible plain- 
text, the procedure is repeated using further trigraphs, but this would oceur only 
very rarely. 

In an example containing 100 cipher groups (m = 100), 50 trigraphs (k = 50) 
were tested, giving some 10,000 trials to be examined. Of these, 280 gave apparent 
settings with 320 (a b c) possibilities. There were 20 causal settings. 

In general, the longer the cipher-text the fewer thé number of trigraphs that 
need be tested. Discarding cases where D or D’ is 0 or 13 mod 26 causes no difficulty 
since sufficient trigraphs are used. The time for testing the 10,000 trials was ap- 
proximately 30 minutes. 


North Carolina State College 
Raleigh, North Carolina 
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Stability Analysis and Integration of the Viscous 
Equations of Motion 


By L. Filler & H. F. Ludloff 


1. Introduction. Techniques are established for the numerical solution of the 
time-dependent, one-dimensional equations of motion of a viscous, heat-conducting 
fluid. The equations of motion are approximated by suitable finite difference equa- 
tions and the stability of the difference equations is investigated by using von 
Neumann’s error analysis. Both explicit and implicit finite difference schemes are 
studied. The implicit equations are solved by an iteration scheme that is formulated 
with the requisite that its convergence does not place conditions on the mesh width 

atio 
ratio 7 

The dimensionless equations of motion to be treated are written in divergence 

form [1] 








dp , OM 

— = 0 
ot Ox 

OM oR _ 
ot Ox 
am, ar _ 
ae 


where p is the density, 1 the mass flow vector, E the total energy, and R and T 
are vectors which are certain non-linear combinations of p, M, and E, namely: 





3—yM _4_ aM/p) 
er ee er ~ 


ME y-1M (; 7) M@M/p) vy AE/p) 
T=y —(5-— -)u————— - - 1 -— 


p a p" . @ p Ox o Ox 


R=(y-DE+ 














with the usual notation that y is the ratio of specific heats (constant = 7), u the 
coefficient of viscosity which is assumed to vary linearly with temperature, and ¢ 
the Prandtl number (constant = ?). 

This system of differential equations is replaced by a system of difference equa- 


sid . _ , 
tions as follows. Derivatives with respect to time 2 are approximated by the for- 
ward difference quotient 
1 
At { fm.n+1 Sm.n} 


where fm,» denotes the value of f at the lattice point z = mAz and t = nAt. Space 
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Listed 6 le 
derivatives IS 


at are approximated by the difference quotients 
rv 


k . 
(Ax)~ > (—1)* (*) Sk’ r,n+e 


=0 


r 
index w = 0 for the explicit difference scheme formulation and w = 1 for the im- 
plicit scheme. Note that this approximation yields backward first differences in 
space and centered second differences in space. 

Both the explicit and implicit difference approximations are used to compute the 
formation of a shock wave from a finite amplitude compression pulse. The computa- 
tions were carried out on an IBM 650 digital computer at the Watson Scientific 
Computing Laboratory of Columbia University. 


k , , ; : , k 
where ( ) denotes the binominal coefficient and k’ is the largest integer < 5° The 


2. Explicit Difference Scheme. The explicit difference equations are with the 
notation Afm—in = fan — fm—-.n 


At 
Pm.nti = Pm.n —~ Ax AM n-1.0 
At 
i s.s ™ Mn . Se ARn- nm 
7 Ax . 
At 
En n _ En eS fo re n 
n+l ° At 1, 


where, for example, 


AR. s.0 — (y i ¥ 1L)AE 1,7 + : ooh A (“) ; 
m—1,n 


9 p 
1 A 
)..j-tmee()., 
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4 M 
— — |Aum-i.n| | A 
3 |Au vl (3 


The linear variation of the coefficient of viscosity with temperature is written 
in terms of the dependent variables as 


tan = (7 — 1) (2) " yaey ] 
p mn 2 p m,n 


The computation of the explicit difference equations proceeds directly, since all terms 
on the right-hand side of the equations are known values taken at the previous time 
cycle. The unknowns on the left-hand side are computed at each lattice point of the 
mesh in the new time cycle. 

The analysis of the stability of the difference equations follows the procedure of 
von Neumann [2]. However, the original system of equations is non-linear, and the 
variational equation obtained for the error has non-constant coefficients. In order 
to apply von Neumann’s method of stability analysis the coefficients are sys- 
tematically set constant, and all derivatives that appear in the coefficients are set to 
zero. Setting the derivatives that appear in the coefficients of the variational equa- 
tion to zero is arbitrary. One could, with as much justification, consider them con- 
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stant. However, the stability analysis is greatly simplified if only those terms not 
having derivatives as coefficients are retained. This does not seem to impair the re- 
sults of the analysis when checked against actual computation. Indeed, when such 
a treatment is applied to the non-viscous hyperbolic equations of motion, the cor- 
rect result is obtained: that stability is determined in terms of the characteristic 
directions [3]. Thus, the variational equation obtained from the system of difference 
equations may be written 


| aon = {I a AlU n,n + BU asin + CU a-in 


where 

p 

Un. =| M 

E | =.9 
Substituting for the error the Fourier term ee*‘e“” leads to the result: 
(2.1) AU, = [J — A(1 — cos 8) + iD sin Uo. 
Here =e and @ = BAz 
and the matrix D = B — C. Note that A = B+ C. Let H = —A(1 — cos@) + 


iD sin @ and write the spectral equation (2.1) as 
[H — KI|JU, = 0 
with K=,i-1. 


: A as. , : : 
The elements of H are, with T = A (1 — cos 0), y = Ar 6, and in combination 
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Here v is the kinematic viscosity coefficient u/p. Now, K the eigenvalues of H satisfy 
the cubie equation 


— [he + hs3|K* — [hyha + heshe — heohss|K — [hi2hesha — hyhayhss} = 0 
and it may be verified by direct substitution that 


‘ ; 1 .~ 2 | 
* = — } = _ = —_— 
(2 2) K ¥ 133 [we ae aha A 


is a root of the cubic equation if ¢ = 3. Factoring this root, the roots of the remain- 
ing quadratic are 





Ee 
a K = —| u@ Par [ 252 yY¥v\ 
(2.3) E +t rls 4/ce +(*2)1 
where 
a = 7(y —- »[F- |. 
Thus, the requirement for stability |\ | = | K + 1| < 1 yields respectively, 
(2.4 \1—|udt2—rii<i 
(2.4) | E + 7 al =< 
| — ie. ae 9 
(2.5) j}1— {ue +X r+ 4/ / 2 +(*2) 2s/<1 
\ ov o Ax . i 
and 


: Y» 1») pel 
(2.6) 1- {ue +-—T- ae 2 hl 51 = 
| ote ~ e+ ms 
Each of the above inequalities must be considered to determine the conditions 
for the numerical stability of the difference scheme. For example, it may be shown 
that inequality (2.4) has a simple geometrical interpretation, viz., the absolute value 
of the locus of the ellipse in the complex 1, é plane is less than or equal to one, where 


A=nt+t 


. i 
n=1-riut22] 


and — = uy. 
Thus, the stability condition obtained from (2.4) is 
At 1 
— 
(2.7) Az ~ 9 : 
utes 


The form of the inequality (2.5) suggests the conservative approximation 
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which gives the stability condition 
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The inequality 2.6 is approximated by 








1 
, 3 a 
Az 
This condition is critical when u = 0, since it can be fulfilled only if 
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Im 41 += >—0. 
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Consequently, we require 
yY v 
cae >a | so . 


For practical purposes, it has been verified when calculations were carried out 
that a reasonable estimate is 


(2.9) + 


ap” ae. ee 
(2.8) Az = 


which together with (2.9) give the conditions for stability of the explicit difference 
equations. 

It may be noted from (2.2) and (2.3) that in the limit as the viscosity tends to 
zero the eigenvalues K degenerate properly so that the correct stability conditions 
for the hyperbolic equations of motion are obtained. 


3. Implicit Difference Scheme. The implicit difference scheme equations are 
written out in full 
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Note that the implicit difference equations are written with the assumption of 
constant coefficient of viscosity u to facilitate solution by an iterative method. 

An inspection of the implicit difference equations reveals that the continuity 
equation occupies a preferred position in the system of equations. The pm,.+: to be 
computed depends only on unknown M’s at the lattice points m,n + 1 and m — 1, 
n + 1. The computations are started from a left-hand boundary, or from a region 
where conditions are known to be constant, so that the point m — 1, n + 1 falls 
on the boundary or in the region of constant value. Thus the values of p, M, and E 
at the lattice point m — 1, n + 1 are known. Then, fer the continuity equation, 
the only unknown is the M,,,.4:. To obtain an approximate value of M4: the 
momentum equation is set up for solution by an iteration procedure. When the 
M m+ has been approximately determined, it is inserted into the continuity equa- 
tion and an approximate value of pm»; may be directly computed. Proceeding in 
this manner along a line of computation m + 1, the Mn4i..4: is next calculated by 
iteration, after inserting the previously determined values of M and p at the lattice 
point m, n + 1. In this way approximate values of M and p are computed at all 
lattice points on the line n + 1. When these values have been determined, the energy 
equation may be set up as a three-point recursion formula for En—i..41, Em.n4i1, and 
Em4i.n41 Which may be solved by a standard reduction procedure. When the first 
set of computations has been completed we have approximate values of p, M, and 
E at all lattice points on the computation line n + 1. The calculations are repeated 
until convergence is achieved.* 

The procedure for computing the momentum equation takes the form of a sys- 
tematic relaxation method. Setting up the difference operator 


wi 4. me x) 
Cows = sage ai — 
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* This procedure was adopted upon the suggestion of H. Keller, Institute of Mathematical 
Sciences, New York University. 
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~ { (EL wt — M...) 


At 
Pe ye i,n+1 


3- M \ 
+2578 ,(F) 
2 | ( Pp | 
where the iteration number is given by j. We then write the following equation for 
Mana 
(3.1) Manet = Mansi + 8Gnnss - 


Approximations to the unknowns at iteration number j are inserted into the ex- 
pression Gi.4: leaving a remainder or residue, since the exact values satisfy 


Ginn+i = 0. The residue is systematically reduced by the factor 6 until convergence 
is obtained. 


The remaining equations are written so that they may be solved directly once 
the approximations to the momentum equation are available. 
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In (3.2) the M2), .4: is known and the M2"! ,, is inserted from the solution of 
(3.1). In (3.3) the coefficients of the ee may be computed after the values of 
M and p have been found from (3.1) and (3.2). The right-hand side of the equa- 
tion is similarly given. Equation (3.3) is thus a linear system of equations in the 
unknowns E2*}..4:. This system is tridiagonal and may be written as follows: 
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The matrix (a,,) is triangularized by letting 
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and 
, 
oy @ hy =e; p #1. 
@(p—1) (p—1) 

Equation (3.4) becomes 

k 
(3.5) Las Bent = be; p= 1,2,---k 

= 


with the triangular matrix a,,-. Equation (3.5) is easily inverted to yield the 
en . 


The convergence of the iteration procedure used for the momentum equation, 
(3.1) Maha = Mans + Gann 


follows the techniques used to analyze the stability of the explicit difference scheme. 
Equation (3.1) is “linearized” and a von Neumann analysis of the “error”’ is made 
[4]. 

Convergence of (3.1) implies convergence of the continuity and energy equa- 
tions. To establish convergence of the iteration scheme we consider convergence 
of the M,,..4: iterates with predetermined values of p and E obtained at each lat- 
tice point from (3.2) and (3.5) respectively. 
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; + At 
ME = Mina + {8 = _ OP ba ss 
ual : Uu <= = A’ Pu-in+i = c= - Ma.n) 
ae 
- (y — 1) ia AE*-1.n41 
(3.6) ; 
At ; 
+ (3 = 7) Uu a AMi-1.041 
3—7 2A ; 
= 9 u Ar dohsnss |} 
MV 
Here vy = =, u = — 
p p 


Before introducing an error for the iterations we note that M,,,. is considered an 
errorless input, and by virtue of holding p and E£ fixed at each lattice point they will 
also be errorless inputs. Therefore, the error term considered is 


Mi, nas - e%e**, 
This term is introduced into (3.6) with fixed p and E, giving 


e* = p + qceos@ + irsin 0 
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8 » At 
[$2 +0-nu] 


._ Af 
6(3— yu. 
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Since q > r° the convergence requirement is 
(3.7) lp|+\q| 
which yields upon substitution of p and q 
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if u = 0. The inequality (3.8) identically satisfies (3.7). However, (3.7) may be 
satisfied with 


1 
(3.9) oe see 
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If (3.9) holds we get from (3.7) that at most 
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1 
1. Tier At 
s+[82+@-u]% 
The least restrictive condition on 6 that covers all cases is obtained from the in- 
equality (3.10). 
Thus, the mesh width ratio may be arbitrarily chosen for the implicit scheme 


but convergence of the iteration procedure is obtained by choosing 6 according to 
(3.10). 


lA 


(3.10) 


4. The Formation of a Shock Wave. The formation of a shock wave from a finite 
amplitude pulse has been numerically treated in three ways. First, a technique due 
to Lax [5] was used, where the equations of motion of an ideal fluid are used to com- 
pute flow fields in which shocks may develop. Second, the explicit difference formula- 
tion of the Navier-Stokes equations has been applied to the problem, with restric- 
tions on the mesh-width ratio as previously determined in Section 2. Third, the 
implicit difference scheme approximation with the previously described iteration 
technique has been applied. An isentropic initial field is prescribed consisting of two 
homogeneous states of different velocity, pressure, and density connected by a 
simple compression wave. The problem is then formulated by asking for the develop- 
ment of this field in time as governed by the equations of motion. In particular we 
look for the time required for the shock to form, the shock’s final shape, and espe- 
cially how the entropy profile, characteristic of shocks, develops. 

The numerical values used are 


(2) — a= 1 p2 = 3.35 
»* M,=0 M, = 452 
7a & E, = 1885 E, = 126 
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Fig. 1.—Non-Viscous Case. 
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Fic. 2.—Entropy Growth Non-Viscous Case. 
The initial velocity decrease from state (2) to state (1) is linear. The unit 
length, ¢, was subdivided into 32 space intervals Ax. For the viscous equations a 
. . v . . . . . . . . 
viscosity parameter -* 1 was taken for both the explicit and implicit difference 
x 


schemes. This value is seen to satisfy the inequality (2.9) with y = 4,0 = ? and 
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Fic. 4.—Entropy Change Viscous Case. 
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Fic. 6.—Entropy Change Viscous Case. 
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a|.-0 = 1. The stability criterion (2.8) for the explicit scheme with the specified 
initial field values is 


= 0.155 for the viscous equations 


RIE RIE 


= 0.366 for the ideal equations. 


For the implicit difference scheme we obtain* 


6 = 0.167 


from (3.10) and arbitrarily chosen x = 1. 


The computation of the viscous implicit equations is somewhat arbitrary; the 
machine controls should be flexible so that the most economical way of solution 
may be found. The values of M and p should first be computed at each lattice point, 
then inserted back into the equations for M and p and computed a second time. 
Then these second computed values of M and p are used to compute EF. The process 
is repeated until convergence is reached. Thus, two “inner” iterations for M and p 
per “outer” iteration for E was found to give the most rapid convergence (9 outer 
iterations per time cycle). The results of the machine computations are shown in 
Figures 1 and 2 for the ideal equations, Figures 3 and 4 for the viscous explicit 
equations, and Figures 5 and 6 for the viscous implicit equations. In each case, the 
time development has been continued until the form of the profiles becomes sta- 
tionary, i.e., only a translation takes place. This was coincident with the time at 
which the entropy had built up to the value to be expected from stationary shock 
theory for a shock of the given pressure discontinuity. 

The velocity profile, which falls off linearly within 32 mesh intervals in the 
initial field, steepens rapidly, especially in the ideal rather than in the viscous case. 

In the ideal case it takes approximately 100 time cycles, in the viscous explicit 
case 140 time cycles, while in the viscous implicit case approximately 21 time cycles 
to reach the final state. But the absolute time for the shock formation, while essen- 
tially equal in the viscous cases (tf = 0.67 explicit; ¢ = 0.66, implicit) is longer in 
the ideal case (t = 1.14). This is, of course, to be expected since the entropy buildup 
is proportional to the viscosity present. 

Note also that the slope of the entropy curves is steep in front, as it should be, 
but falls off to zero in the rear, indicating that the air particles at the rear have not 
been swept over by the fully developed shock. 
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*§ = 0.2 was found in practice, to give convergence of the iteration scheme. 
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A Simple Experimental Computer with 
Negative Basis 


By A. Lazarkiewicz and W. Balasifiski 


Summary. This paper presents the technical data, logic, and control organiza- 
tion of a simple experimental computer operating in the minus-two system. The 
principles for composing the instructions from elementary operations, the char- 
acteristics of minus-two computer arithmetic, and the logic of the arithmetic unit 
are briefly explained. The paper is divided into the following sections: 1, Introduc- 
tion; 2, Fundamental Logical Circuits; 3, The Memory; 4, Block Diagram; 5, The 
Control Unit; 6, Coordination with Teleprinter; 7, Instructions; 8, Minus-Two 
System and Range of Numbers in the Computer; 9, The Adding Unit; 10, The Sign 
and Overflow Register; 11, Multiplication; 12, Acknowledgments. 


1. Introduction. An experimental negative-base digital computer is being tested 
at the Warsaw Technical University. It is a small general-purpose machine for sci- 
entific and technical calculations which works at a speed of about 100 oper/sec. The 
computer makes use of the serial mode of operation; it is a single-address, fixed- 
point machine with a magnetic drum as the only memory. Punched paper tape has 
been used for input of data and the results are printed by means of a standard 
teleprinter. 

It is intended to supply some university computation centers in Poland with an 
improved version of this computer for training purposes. 


2. Fundamental Logical Circuits. Fundamental logical circuits have been 
worked out in dynamic tube technique with the use of germanium diodes for obtain- 
ing “and” and “or” logic. The arithmetic part of the machine and its control have 
been built up from a few normalized types of these circuits. The circuits are shown 
in Figure 1. 

The delay unit is based on the Havens circuit used with different modifications 
in computers like the Pegasus, NORC, in the IBM 701 to some extent, etc. After 
an input pulse has been applied, an impulse, regenerated both in time and ampli- 
tude and delayed by a clock unit of time, appears at the output of the delay unit. 
A dynamic trigger may be obtained by feeding back the output to the input, as 
shown in Figure 2. 

At the output of every fundamental circuit there is a cathode follower. A pulse 
representing the digit ‘‘one’’ has about 30 volts amplitude and lasts approximately 
7.5 ws. The fundamental circuits are linked into blocks according to fixed rules 


which take into account delays and pulse deformations in amplitude caused by gates 
and negations. 


3. The Memory. A magnetic drum serves as the memory. It is made of alu- 
minium covered with iron oxide about 20 » thick. The drum is fitted with heads set 
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Fic. 2.—A Dynamic Trigger. 


about 30 » from its surface. They are made of permalloy strips 1 mm wide on which 
160 turns of wire are wound. Over 2700 rpm give a maximum access time of 22 ms. 
A read-record technique of phase modulation has been adopted. A “one’’ is recorded 
by a current pulse flowing through the coil of the head, immediately followed by an 
equal but oppositely directed pulse. The direction of both current pulses is reversed 
in “zero” record. Each zero-one sequence contains an even number of positive and 
negative pulses so it has no direct component. The read-frequency varies from 50 ke 
for 010101 --- to 100 ke for ‘‘zeros”’ or ‘‘ones”’ alone. 

Switching the tracks for reading and writing is accomplished by means of an 
electronic switch. The organization of the computer allows neglect of the transition 
states in read control circuits. 

In addition to the memory itself there are 3 dynamic circulating registers on 
the drum, each about 1 word long. They are the accumulator, the current order 
address register (the so-called counter register), and read-out instruction register. 

Drum registers are equipped with a writing head and a reading head which are 
placed close to each other on the same track, as well as amplifiers for writing and 
reading which contain pulse-forming circuits. The two heads can be positioned to 
secure correct performance. To eliminate unwanted signals which may be picked 
up, the heads of the registers are screened with permalloy strips placed inside the 
cylindrical case holding the head (the case diameter is 7 mm). 

Six-bit addresses of separate words are written on the address track and suc- 
cessive addresses are placed every 13 words, as below 
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Fic. 3.—Block Diagram of the Computer. 
Symbols: 
z—drum memory a—accumulator register 
r—instruction register d—input-output register 
l—counter register >-—adding unit 
v—track register w—sign-overflow register 
s—control register Ro , lo , Xz, Az , et cetera—gates 


This system of addressing has been chosen in order to facilitate optimum program- 
ming work in respect to the access time. 
The capacity of the memory is 4096 words. 


4. Block Diagram. The block scheme of the computer has been so designed that 
it is possible to combine transfer of numbers and instructions with performance of 
operations. 

The adding unit is the central element of the transfer (Fig. 3). It is linked with 
the counter register “J,” the instruction register “r,’’ the memory and the input- 
output register “d’”’ on one side, and with the accumulator “a” on the other side. 

Writing on the drum is done directly from the accumulator or from the counter 
register. Individual parts of the computer are linked through control gates. To 
execute an instruction, the appropriate gates on the transfer line should be open, 
e.g., the line of transfer through gates Xo and Az must be free to bring a number 
from the memory into the accumulator. To add “one” to the counter register “1,” 
gates Lo, Lz, Ns must be open. Gate Ns makes it possible to introduce plus one into 
the adder as the third argument in any operation. 

The accumulator register is a drum register of one word length with the possi- 
bility of right and left shift. Figure 4 gives a simplified scheme of this register. 

The basic loop is closed through two delay-lines, and the drum delay T — 2 
equals 34 clock periods. The introduction of a new word from the adder takes place 
at Az = 1, and it automatically erases the previous value by breaking the loop by 
means of the negation n . Shifting to the right occurs at Ap = 1; it shortens the 
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p 4 drum delay T-2 . 























to= 
Fie. 4.—Logical Circuit of the Accumulator Register. 


loop of the accumulator for a word period by one bit; the normal path is then 
blocked by negation n; . When shifting to the left (Al = 1) the loop is lengthened 
by one bit; its normal path is blocked too. 

The “‘d” register is a 5-bit trigger register serving as a buffer between the tele- 
printer and the computer. Its five positions correspond to five holes in the punched 
paper tape. Loading and unloading the “‘d’’ register by a converting device coor- 
dinated with the standard teleprinter is accomplished by means of a matrix con- 
trolled by a separate counter. 

The current address of an actually executed instruction is kept by the counter 
register; the contents of the counter corresponding to the operational part of the 
order equal zero. During multiplication the counter register is used as the multiplier 
register. 

The ‘‘r” register is the instruction register. It serves to keep order during the 
waiting time required for access to the right memory address from which the num- 
ber is to be read out, and it acts as a buffer for the control register. During multi- 
plication it stores the multiplicand. 

The “v” register is a six-bit trigger register which holds the number of the track; 
it controls the matrix which defines the proper track on the drum memory. 

The control register “‘s’’ holds the operational part of an instruction. It consists 
of twenty triggers, the states of which are determined by a serially introduced com- 
mand from the ‘“‘r” register or by parallel setting from the simple matrix of the con- 
trol unit. 

The “one” state of each position of the ‘‘s” register directly opens a gate for 
each appropriate elementary operation. Positions Xo and Rz are an exception, and 
their active states break their corresponding paths while “‘zero”’ states open them. 


foput —> |Xo] Do|Adt|Ro| Lo | Ne| Qe 


| parallel input { 





Ao 








At \A> 





Ne| Sg 





We| Mo} Xz| De| Ae|L|Oe| A 











Fie. 5.—The Control Register. 
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LUT TT TTT TT 


+ 9perational part ——— 20 ———e{ |e 1th 6 -o\-— position 6 -» 
k— address part —12 —» 


Fic. 6.—The Structure of an Instruction. 




















5. The Control Unit. An instruction is defined by 32 bits placed as follows (Fig. 
6). 
Twelve digits are reserved for the address part, 6 bits point out the proper 


track, and 6 bits—the time (position on the track). The operational part has 20 
digits. 





The machine works in three steps; namely, 1) forming the next instruction 
address, 2) extracting the instruction from the memory, 3) the execution of the 
instruction. 

1. Instructions are placed in the memory in successive addresses. Defining the 
new instruction address requires the addition of “one” to the counter register “1.” 
This is done by the adding unit. The new instruction address is forwarded at the 
same time to the “r” register. That occupies a period of time of one word, i.e., 36 
clock units of time. During the next word period its operational part is shifted to the 
““s” register and the address part, defining the track on the drum store, to the “v” 
register. Since there is no “1” in the contents of the counter (in its operational part) 
we shall find only zeros in the “‘s” register, which means the operation of the second 
step. The time part (position part) of an instruction circulating among the others 
in the “r’” register is compared to the output of the address track. As soon as 
coincidence is reached the second step begins. 

2. In the second step the instruction is read out from the memory and brought 
through the adding unit to the instruction register ‘“‘r.” Then the operational part 
is shifted to the control register ‘“‘s,” and track address to the track register “‘v” in 
the next word period. Similarly, the time part of the address circulates in the “‘r” 
register and is compared with the address track. Coincidence starts the third step. 
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Fic. 7.—Block Diagram of the Control Unit. 
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3. During the third step the operation defined by the contents of the “‘s’’ register 
is executed. 

It is evident that the minimum execution time for one instruction is equal to 
five word periods, i.e., 1.8 ms. 

If the instruction to be obeyed is conditional and the sign-overflow register is in 
the “‘one”’ state, the “end of cycle’’ pulse is produced just before the execution time 
(the third step), which clears the control register and places the elementary opera- 
tion of the first step in it. Thus, each (conditional) instruction can be ignored, de- 
pending upon fulfillment of a given condition. 

The machine may be stopped by a stop-order or during reading and printing. 
The computer is also stopped in the cycle-by-cycle mode of operation before the 
execution of the instruction placed in the control register ‘‘s.”” Step-by-step opera- 
tion causes a stop after every step. 

To start the program from the address n, the number m must be placed in the 
counter register “1,” e.g., from the ‘“d’’ register by means of manually operated 
switches. 


6. Input-Output. The “d’’ register and its converting unit link the computer with 
the input-output units (the teleprinter and the reader). Five positions of the “d’’ 
register correspond to the holes in the punched paper tape. Writing or reading is per- 
formed by the converting unit in a parallel manner with a matrix controlled by a 
separate simple counter. The converting unit produces the end-of-reading and the 
end-of-printing pulses according to the following rules, and receives the start of the 
printing pulse. 

1. The instruction ‘“‘read from d” is stopped until the end-of-reading pulse comes. 

2. The instruction “put in d’’ sends the print pulse to the converting unit. Dur- 
ing printing the next “put in d” instruction is withheld. 


7. Instructions. The instructions consist of 20 elementary operations defined 
by the positions of the control register ‘“‘s.’”” The function of these positions is directly 
visible from the block diagram. Several remarks should be added: 

1. Adr—the operation to which this position is adjoined is executed only on the 
address part of the word. 

2. Mn—position “multiply” initiates the cycle of operations forming the 
product. 

3. In order to write the contents of register “7’”’ in the memory the positions Ro 
and Lo must be activated. 

4. The “stop” order is realized by the positions Do = Dz = 1. 

5. The execution of an instruction to which position Wa is adjoined depends 
upon the state of the sign-overflow register. If the register is in the “1” state, this 
order is omitted and the computer performs the next instruction. 

The number of instructions composed of elementary operations is very great. 
Several operations may be executed simultaneously, and this increases the effective 
speed of computing. 

Examples of instructions: 
(x + aja n—add the number from address ‘“‘n” to the contents of the 
accumulator 
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aL —shift the contents of the accumulator one place left 

Wa Adr (r)jl n-1—if the sign-overflow register is in “0” state put the address 
part of this order into counter register “/’’ (jump to the 
instruction in “n’’); if the sign-overflow register is in the 
“1” state carry out the next instruction 

Adr Sg (Nsr)a n—add 1 to the address part of the instruction in the “r’’ 
register, transfer the result into the accumulator and 
examine the sign 

Adr (1)z, (r)l n—put the current instruction address in memory cell “n’’; 
set ‘‘n’’ in the counter register (jump to a subroutine) 

aL, (Nsr)r (10924-n)—shift the contents of accumulator “‘n” places left (n = 2); 
instruction aL is repeatedly executed until the carry pro- 
duced by the address part of the instruction blocks the 
position Rz in the register “s.” 

A simple method for placing the input routine in the memory as well as a system 
for coding instructions has been worked out. A special routine located on a few pro- 
tected drum tracks makes it possible to prepare programs even without knowledge 
of the internal organization and the number system used in the computer. 


8. Minus-Two System and the Range of Numbers in the Computer. Each real 
number may be represented in the minus-two expansion 


N 
@ = (Gq +++ Gy0o-_40_2 -*+)-2 = >, a;-(—2)' 
3 


where a takes the values of 0 and —1 written as 0 and 1 (the so-called negative 
notation) and N is the greatest i for which “‘a;,”’ still has a nonzero value. 
In the minus-two system the weights of expansion orders are as follows 





expansion order (i) --- 3 3.1 0 -1 -—2 -3-:-- 
value of ‘1” 8 -4 2 -l , —1 2 


For expansions of positive numbers N is odd, for negative numbers N is even. 
For numbers limited to imax (imax 2 N), the position of zero is nonsymmetrical, 
e.g., When imax = 1 the numbers are comprised within the range 


1 
—$ Sa}. 


Shifting a number one place left or right is equivalent to multiplication or di- 
vision by 2 with inversion of sign. 
Examples of numbers represented in the (—2) system: 


—3 111 —4 0.010101 --- or 1.101010 --- 
—2 110 —- 0.010000 --- 
—1 1 — 0.011111 --- 
0 0 0 0.000000 --- 
1 11 t 0.110101 --- or 0.001010 --- 
2 10 i 0.110000 --- 
3 1101 4 0.111111 --- 
4 1100 3 0.100000 --- 
5 1111 3 0.101010 --- or 11.010101 --- 
6 1110 3 11.010000 --- 
7 1011 1 11.000000 --- 
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The algorithms of operations a + 8B, a — 6 and —a + 8 include two carries 
~ = +1 and p, = —1, and the algorithm of operation —a — 8 only one. The adder 
for carrying out this algorithm has a structure equivalent to the adder of the kind 
a + 6 for the well known binary system (+2) (see Table 9b). 

Multiplication can be performed by successive subtractions of the multiplicand 
from the shifted partial product according to the values of digits of the multiplier. 

Algorithms for other operations such as division, rounding off, symmetrization 
of the range and even finding square-roots have structures similar to the algorithms 
in the (+2) system commonly used and require an almost equal quantity of 
hardware. 

Due to the experimental character of this machine, only multiplication was ac- 
complished; other operations are left to programming. 

The 36-bit word adopted in the computer is interpreted as follows 


@j00°A1A_2 **- A_330_-34. 


In the memory, however, only digits behind the point are held, i.e., the numbers 
in the range 


—t <a < 3. 
The error caused by rejecting further bits lies in the range 
—0.19-10°" < 6 < 0.39-10™. 


The accumulator contains the whole word, e.g., four times the greatest value of 
a number read out from memory. 
In multiplication both arguments can be represented as 


Odo-@_10_2--*- Gy. - 


This enables execution of many operations without overflow. 


9. The Adding Unit. The adding unit is composed of two sign-inversion circuits, 
P, and P2, as well as the adder proper. 

Units of sign inversion are controlled by positions 0a and Od of the ‘‘s’’ register. 
When 0a = 0d = 0 (i-e., 0d = Oa = 1) both circuits P invert the signs of argu- 
ments a and @ and the sum y = a + 8 appears in the output of the adding unit. 
If 0a = O and Od = 1 the sign of the number 8 is not changed and y = a — 8. 





Fic. 8.—Block Diagram of the Adder. 
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a; Pi (—@); Ris 
0 0 0 0 
1 0 1 1 
0 1 1 0 
1 1 0 0 





Fic. 9a.—Table of Sign Inversion of the Number a in —a. 














' a; b Pi Ci Pin | 
0 0 0 0 o | 

1 0 0 1 1 | 

| 0 1 0 1 1 | 

1 1 0 0 1 

0 0 1 1 0 

1 0 1 0 0 

| 0 1 1 0 0 

1 1 1 1 1 








Fic. 9b.—Table of Operation y = —a —8. 
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Fie. 10.—The Circuit for Inverting the Sign of a Number. 




















Position Ns of the control register adds one to the result of normal addition, 
eg.,y =a+t+ B+ (—2)™. 

The algorithms of sign inversion and of operation y = —a — 6 are presented 
in Tables 9a and 9b. In these tables a; , b; denote digits of arguments a, 8; (—a);, 
c; — digits of the result —a, y; pi — carries. 

Inversion of the sign of a number is accomplished by the logical circuit shown 
in Figure 10. 

The suppression of carries by the position 0a = 1 has the effect of transporting 
the incoming number without changing the sign, which causes subtraction. The 
existence of an initial carry means adding of (—2)~“ to both arguments but, of 
course, only when 0a = 1 or Od = 1. 

It is evident from the truth-table for the adding unit that if the carries are 
suppressed and p; = a; + b; is assumed (in Boolean algebra), the adding unit per- 
forms the operation of the logical product y = a & 8. In this case the circuits P 
do not invert the signs of the arguments, i.e., 0a = Od = 1. 
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Fic. 11.—The Sign and Overflow Register. 


10. The Sign and Overflow Register. When position Sg in the “‘s’’ register is 
“1,” each number coming out from the adding unit is advanced to the sign defining 
unit. If the examined number is positive or equal to zero the register is left in 
state ‘‘1’’; if not, the register is cleared. For this purpose the clock pulse (marked 
© in Fig. 11) turns on the trigger. Bits with value “1” in even expansion orders of 
the number switch off the trigger (which was previously turned on by the clock 
pulse @)), and in odd orders switch it on. Thus the state of the register is determined 
by the “one” bit of the highest expansion order (including ao). 

The overflow of a number in the accumulator or ef a number forwarded there 
is examined in a similar way. The trigger is turned on if bit a) (the © pulse in Fig. 
11) is equal to “1,” and is switched off in the opposite case. 

The contents of the sign and overflow register are not changed until a new in- 
struction arrives. 

11. Multiplication. Multiplication consists of the preparatory steps, proper 
multiplication, and final steps. In the last clock time of every step new positions 
defining the next step are set in the “‘s” register with a simple matrix. This matrix 
is controlled by the special positions of the “‘s” register and by the coincidence of 
suitable addresses taken from the drum track. 

One of the arguments of multiplication (the multiplier) is placed in the ac- 
cumulator, the other (the multiplicand) in the ‘“‘r” register; when squaring the 
argument is placed simultaneously in “r” and “a.” 

Counting the multiplication steps begins with address zero of the address track, 
and the contents of the counter register “l’’ are transferred to this location (track 
zero, position zero). 

Simultaneously the multiplier is placed in register “l’’ and the accumulator is 
cleared. Then proper multiplication follows: the contents of ‘‘l” are shifted right, 
the outgoing “1” bit of the multiplier opens the path for the multiplicand (with 
inverted sign) to the adding unit, the contents of “‘a” are shifted right, the outgoing 
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bit is placed in the emptied position of the register “1.” Coincidence of the constant 
address with the address track ends this procedure. 

After completing multiplication the contents of register “l’’ (the least signifi- 
cant part of the product) are written in the storage cell on track zero and in actual 
position under the writing head. Finally, after the waiting time is completed the 
contents of the zero address are brought again to register “1.”” Consequently, the 
most significant part of the product is to be found in the accumulator, and the 
least significant part in memory cell —7. 

The speed of multiplication is 16-36 mult./sec., depending on the access time. 

This scheme makes it possible to perform many operations with a single in- 


: 2 a 
struction, e.g., x-2, ©-a, a-a, —a-a, 4a, ru (x + a)-a, (x — a)-a, (x + a)’ ete. 
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TECHNICAL NOTES AND SHORT PAPERS 


The Approximate Solution of an Integral 
Equation Using High-Order Gaussian 
Quadrature Formulas 


By Stephen M. Robinson and A. H. Stroud 


The purpose of this note is to illustrate the usefulness of the high-order Gaussian 
quadrature formulas given in [1], [2], [3] in a problem of approximating the value of 


(1) I =| g(x)oe(x) dx 


where ¢(2x) is the solution of a Fredholm integral equation 
b 
(2) o(z) — [ K(2, wely) dy = f(z). 


Here K(x, y), f(x) and g(x) are known functions and K(z, y) has a discontinuity 
at y = 2. 

The approximate solution of (2) was found by a method described in [4], p. 101. 
It consists of replacing (2) by 


ota) [1 ~ [ K(a,v) ay | — | K(z, y)loly) — oa) dy = f(z). 


b 
Now, since K(2, y) is a known function, / K(x, y) dy can be determined as a 
a 


function of xz, say x(x). It is reasonable to suppose that the integral 
b 
(3) [ Kz, Whey) — o(2)] ay 


can be better approximated by a quadrature formula than the original integral in 
(2). If we use a quadrature formula with weights A; and points x ,k = 1,2, --- ,n, 
to evaluate (3), then approximate values of the function ¢(2;,) may be determined 
by solving the linear system 


¢(2x;)[1 — x(x)] — > A, K(2;, te) [e(a%e) — 9(%i)] = f(@i) t= 1,2,---,n. 


The same quadrature formula may then be applied to approximate the integral (1). 

To illustrate the accuracy of Gaussian quadrature formulas in this problem we 
give the results of some calculations concerning flow of a gas through a conical 
opening which were carried out for Prof. R. P. Iezkowski of the University of Wis- 
consin Extension Chemistry Department. The kernel K(x, y) was defined for two 
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TABLE 1 
Values of I 

s a=0 | a= i10 | a= 10 

L=2 L=2 L=10 

“ } a = 
12 0.5142287 _ _ 
16 0.5142294 | _ — 
24 0.5142305 0.6705783 0.4768471 
36 0.5142306 0.6706079 0.4775303 
48 0.5142306 0.6706184 0.4777735 
al 
K6,Y) 





pe — 
o i123 45 6 78 9 10 


Fic. 1.—Cross-Section of the Kernel K (x, y). 


parameters, a and L. The table gives values of J for three cases; in each case J 
was evaluated using the indicated n-point Gaussian quadrature formulas. The 
calculations were carried out using single-precision, floating point arithmetic. 

The nature of the kernel is illustrated in the figure which gives a typical cross- 
section of K(x, y). 

The calculations for a = 0, L = 2 area good illustration of the accuracy achieved 
over most of the range investigated: a = 0 to 89, L = 0.1 to 10.0. In this case the 
24-point formula seems to give 6-decimal accuracy. The other cases of the table 
are some of the least accurate. For a = 10, L = 10 the 24-point formula seems to 
give 3-decimal accuracy and the 48-point formula about 4-decimal accuracy. In 
all of these calculations it appears that, roughly, doubling the number of points 
produces another decimal place of accuracy. No analytical estimates of the error 
were carried out. These calculations illustrate, however, that by using 8 significant 
figures in such calculations it is possible to achieve 6- or 7-figure accuracy in the 
result if a sufficiently accurate formula is used. 

The calculations for n = 12, 16 were carried out on the IBM 650 at the Numeri- 
eal Analysis Laboratory and for n = 24, 36, 48 on the IBM 704 at the Midwestern 
Universities Research Association. The programs were written by Mr. Robinson. 
Funds for the computations were provided by the University of Wisconsin Chemis- 
try Department. 
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Bi-Variate, Rectangular, Optimum-Interval 
Interpolation 


By Ferdinand Freudenstein 


1. Definition of “Optimum-Interval” Interpolation. Let the real function f(z, y), 
defined in the interval —h Sx Sh, -—kK Sy & k,h >0,k > O, and assumed to 
have continuous partial derivatives of order up to (m + n + 2) in the rectangle, 
be approximated by the polynomial 


(1) Prn(x,y) = Dy Ds Grat”y’ 


r=0 s=( 


so that at the (m + 1)(m + 1) precision points, (xu, yo), u = 0,1, 2, --- m; 
9 = 0, 1, 2,---2, 


(2) Plt ; Yr) —= f( Tus, Yr Die 
In that case, 
(3) f(x, y) = Pamn(x, y) + Ran 


where the remainder, Rinn , is given by the expression (see [1]) 





Noga eae an" v) a”*"f(z, 0) 
Ran = Ga papi Lh — %) Geer + G ; Ils I YU — ) 


= grrtoreten ay, ') 
~ weer © -~) Ie - w ogee 


(4) 





where &, &’(, 7’) lie between the greatest and least of the numbers 2, 2,(y, y»)- 
Consider the set, S, of all precision points (2, , y,») within the rectangle (—h, h) 
(—k, k)—called (h, k) from now on. S will be called the spacing. Let h = ké, where 
6 is a positive number defining the proportions of the rectangle. Furthermore, we 
a” f(x,y) a” f(x,y) 


restrict ourselves to functions for which the partial derivatives axe? Gye? 





' Sioa y) 


da™ tl ayntt 
limit on the maximum absolute value of R,,, regardless of spacing. It is assumed 





are not equal to zero within (h, k). Let R be the largest lower 
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that f(z, y) is not of the form P,.»(x, y); thus R > 0. A particular spacing, So, 
will be determined with the property 





i (Mum absolute value of R,.. with spacing S) 
im -= 1, 
h+0 R 

Tu 


6 and all i? ve remaining fixed. The spacing Sp is called an optimum spacing. It 


should be observed that for finite values of h, the maximum absolute value of the 
remainder will in general not be minimized exactly, because of the unpredictable 
nature of the variation of the partial derivatives in equation (4). 


2. Optimum-Interval Determination. Let 


eae. " wee we ee 
." (m+! oer’ "GSD! ot? 














1 "sede tet ane t’n’) 
and Ci,’ = (m + 1)!(n + 1)! OE ™HGy/aH : 
Then 
(5) ct’ Run = P+Q+ PQ 
Cz Cy 
where 





P = P(z, #',n',n) = ee I (x — 2x); 
(6) ‘ 


Q _ Qy, - 1’, é) - oe I] (y — Yo). 


As was assumed above, the partial derivatives are not equal to zero, so that their 
bounds can be identified as follows: 


% S\o|)Sm;hS\eq|) Shs 


lA 
. 
lA 
e 


Since 
(7) II (a — x.) < (2h)"* and JT] (y — y,) < (2k)"™ 
u=0 v=0 
the maximum absolute value, M, of (PQ) within (h, k) regardless of spacing, satis- 
fies the inequality 
2 


(8) M s -2- (2h)""(2k)"". 
a, b; 





It will be shown that an optimum spacing, S, , is that spacing obtained when z, , 1, 
are given by the roots of T,4:(2/h) = 0, Tnis(y/k) = 0, where T;(z) is the poly- 
nomial of degree k in z, with leading coefficient unity, deviating least from zero 
within the interval (—1, 1) [2]. These roots are: 


2a (2u + 1)x. it ae (20+ 1)x 
(9) t= —hoos ray | Y = kes 9 


’ 
t 
? 
; 
; 
3 
2 
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This spacing will also be called Chebyshev spacing. Thus, if W* = | P| + | Q| 


(10) (: omit + 1 ae) < W < (> ia ali + C2 pw) 
by 17) by ay 


v* 


where W is equal to Winax with Chebyshev spacing. 


It follows from relations (8) and (10) that 


= | es 
(11) lim Ww 0, 6 and all x,,/h, y./k remaining fixed. 
h>0 
We can now determine a non-zero lower limit on the maximum absolute value of 
Rinn as follows: 
We restrict ourselves to such rectangles (h, k) for which M/W is less than 
unity, which, using relations (8) and (10), is certainly true when 





1 
aq by C1 ) 1+ where = m if m<n 
B <4 en fe "\ { 
gina \- az + 02.2") | Pen 2 ase 


For any spacing whatsoever, consider equation (5) at the point (2*, y*) at which 
W* has its maximum value within (h, k): 


(12) | Cera’ Fem | = Wmex + \M (0 
| C€gCy | (@*,y*) 


Since the extreme values of | P| and | Q| occur independently, it follows that 
Wrox = Wmin , Where Wmin is the lower bound for W in equation (10), ice., 


y Ci a c —nT n 
W in 2 — + “1 2 ere 


IIA 


r 


lA 


1) 


be Ay : 
Hence, 
(13) PE | Rae ] > Wnin — M > 0. 
| CzCy | (z*,y*) 
We are thus led to define R by the equation 
(14) 2. R = Wain — M. 
ay, by 


We know then that for any spacing, there exists a point (2*, y*) within (h, k) for 
which | Rinnn| 2 R. Let Ro = | Rin |max for Chebyshev spacing. 


M\\ 
“pay [Q) 
(15) lim (%) = lim (y—4) 7 a 


noo \R n0 \Wain — M ho0 (7) ite. 
a“ 


6 and all x,/h, y»/k remaining fixed. Thus, optimum spacing (as defined above) is 
obtained when z,, y» are respectively the zeros of the Chebyshev polynomials 
T m4i(2/h) and Trai(y/k) = 0. 


3. Bi-variate Polynomials P,,,,(x, y) With Leading Coefficient Unity Deviating 
Least From Zero Within (h, k). Let Pmn(x, y), defined as in equation (1), but with 
leading coefficient unity, (an. = 1), be written as follows: 
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(16) Prn(z,y) = 2"bm(y) + 2” “bmaly) +--- + 2bi(y) + do(y) 


where ¢:(y) = >) joay’(k = 0, 1, 2---m) and om(y) has leading term y”. 
Let y = ym* be a value of y for which | ¢,.(y) | reaches its extreme value within 
(h, k). Then 


(17) | Piel &, y) |max > | Pm( Ym*) | | Qual 2, Ynn) jmax 
where 
(18) Qmn(x, Ym*) = 2" + > oe") - oF. 


The least possible values of | dm(Y¥m*) | and | Qmn |max are obtained if 
om(y) = k*T,(y/k) 

Qmn(Z,Ym*) = h”T,,(x/h) 

For this to be the case it is necessary that 


(19) 


1,2 --- (m — 1) 
»2,°°- (n+ 1) 


and a, is the coefficient of x* in the Chebyshev polynomial h”T,,,(x/h). Any of the 
equations (20) is an algebraic equation in y,,* of degree n or less, having at most n 
real roots. However, ¢m(Ym) in equation (19) has an identical extreme value (except 
for sign variation) at least at (n + 1) points. Hence for equation (20) to hold at 
these (n + 1) points it must represent a set of identities, ie. = ak"T,(y/k), 
hence 


> 5 — m 7” x y 
(21) P..(2,y) = h°k'T.. () v. (2). 


This polynomial deviates least from zero within (h, k) with maximum deviation, 
L, given by L = 4h"k"2"‘"*”. Note that the lines x = x, , y = y, have replaced the 
points x, in the one-dimensional theory and that the symmetry P,,,.(+2, +y) = 
Pinn(X,Y) = Pmn(y,2) is preserved. Generalization to n-dimensions follows at once. 

The above developments arose in the course of work on the synthesis of mecha- 
nisms having two degrees of freedom and the author has also utilized with benefit 
and appreciation the works on interpolation by Biermann [3], Hildebrand [4] and 
Kunz [5]. 


$4: (Ym™, ) o k= 0, 
(2) a ™ m* = | 


Columbia University 
New York 27, N. Y. 


1. J. F. Sterrensen, Interpolation, The William and Wilkins Co., Baltimore, Md. 1927 (p. 
206, eq. 13). Reprint available from Chelsea Publishing Co. Formula. goes back to S. Narumi, 
The Tohoku Math. J. v. 18, 1920, p. 313. 

2. N. I. AcHTESER, translated by Cuares J. Hyman, Theory of Approximation, F. Ungar 
Publishing Co., New York, 1956; in particular p. 55-60 on Cheby shev theory. 

3. O. BIeERMANN, Vorlesungen Ueber Mathematische Naeherungsmethoden, F. Vieweg und 
Sohn Verlag, Braunschweig, Germany, 1905; in particular p. 138-146 on interpolation involving 
two “ons 

F. B. HrupEBRAND, Introduction to Numerical Analysis, McGraw-Hill Book Co., New 
¥ ork, 1956; in particular p. 279-282, 386-391 on Chebyshev polynomials. 

'K.S. Kunz, Numerical Analysis, McGraw-Hill Book Co., New York, 1957; in particular, 
Chapter 11, p. 248-274 on interpolation in tables of two or more ‘variables, and equation (11.60) 
on page 269. 
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Divisors of Mersenne Numbers 
10,000 < p < 15,000 
By Sidney Kravitz 


Riesel has previously presented [1] a table of the smallest factors, g, of the 
Mersenne numbers, M, = 2” — 1, p, a prime, subject to the conditions that p < 


Divisors of Mersenne Numbers 
q Divides M, = 2?—1; K = (q—1)/2p 





— 

















K ” | K | t K | » K p K 
10,007 | 12 10,883 | 1 12,037 3 | 12,923 1 | 13,841 3 
10,067 | 25] 10,891; 9] 12,041 | 75) 12,941) 4 13,883 1 ) 
10,091 i 107949 | 7) 12,049 | 11 | 12,959 1 | 13,903 | 12 | 
10,111 5| 10,973| 3] 12,073| 15] 12,979 | 60] 13,963 | 12 
10,133 7| 10,979} 4] 12,107| 4{/13,001| 4/ 13,997; 38 
107139 | 60} 10,987; 9| 12,119| 1] 13,003; 8| 13,999) 5 
10,141 8] 10,993 | 11] 12,157 | 132] 13,049 | 4) 14,033) 3 i 
10,151 | 96] 11,059 21] 12,203 1 | 12,063 | 20) 14,057, 12 : 
10,163 1} 11,071 | 113 | 12,241 | 288} 13,109 | 7] 14,071 | 69 ; 
10,181 | 24] 11,083 | 336 | 12,263 1 | 13,127 | 108 | 14,081 | 15 
10,211 | 4] 11,119! 5] 12,329 | 207 | 13,163 | 28) 14,083 | 5 
10,247 | 9] 11,171 1) 12,347 108 | 13,17 | 3 | 14,149 | 11 4 
10,271 | 1] 11,173 | 12) 12,379 5 | 13,229 | 15 | 14,159 1 
10,289 | 12] 11,197 | 12] 12,391 8 | 13,241 4 | 14,197 | 288 : 
10,301 | 220] 11,243 | 57) 12,401 | 51 | 13,249 | 96 | 14,303 1 I 
10,321 | 3] 11,273 | 27] 12,409 | 396 | 13,291 5 | 14387 13 
10,331 | 1] 11,299] 12) 12,421 3 | 13,313 | 40] 14,411 | 13 
10,333 | 60) 11,311 | 65] 12,437 4 13,339 | 5] 14,461 | 20 1 
10,357 3 11,317 | 108] 12,451 | 200 | 13,367) 4) 14,537) 3 . 
10,429 | 267} 11,399) 60] 12,479| 4) 13,397 | 19 | 14,551 5 : 
10,433 | 7) 11,471 1 12,527 33 | 13,411 | 20 | 14,591 | 25 M 
10,457 | 3] 11,497 | 23] 12,539 | 4 | 13,417 | 12 | 14,5933 
10,513 | 12] 11,519 | 1] 12,547 | 129 | 13,451 | 1 | 14,627 88 
10,589 | 96] 11,579 1 | 12,553 | 32 | 13,463 | 1 | 14,633 | 40 
10,607 | 13] 11,617 | 80] 12,577 | 368 | 13,487 | 4 | 14,639 57 
10,613 | 16] 11,621 | 36] 12,583 | 8 | 13,537 | 80) 14,657 | 3 
10,631 | 25] 11,699 1 | 12,611 | 49 | 13,567 | 377 | 14,669 4 
10,657 | 104 11,701| 3 | 12,641 | 31 13,57 | 3 | 14,699 1 
10,663 | 8] 11,719 | 182 | 12,647 | 340 | 13,619 | 1] 14,713 | 20 
10,687 | 332 | 11,783 1) 12,671 1 | 13,627 | 32] 14,767) 5 h 
10,691 | 1] 11,813) 3) 12,703) 5] 13,669| 11] 14,779) 9 q 
10,733 | 3] 11,827 | 113} 12,739 | 45 | 13,679 | 4 | 14,783 l | 
10,753 | 8] 11,831 1] 12,757| 3 | 13/681 | 44 | 14,831 1 se 
10,799 | 1] 11,833| 3] 12,781 | 48) 13,697| 3 14,843| 16 th 
10,837 | 83] 11,939) 1) 12,791 1 | 13,751 | 9 | 14,879 l ec 
10 ,853 3 | 11,959 | 5 | 12,823 | 20 | 13,759 | 144 | 14,891 | 16 se 
10,859 | 25] 11,969 | 12 | 12,899 1 | 13,763 | 1 | 14,939 l 
10,861 | 68] 11,981 | 171 | 12,911 | 21 | 13,790 | 16] 14,957; 4 
10,867 | 12| 12,011 | 1| | | th 
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10,000 and g < 10,485,760. The purpose of this note is to present an extension of 
this table for 10,000 < p < 15,000 and g < 10,485,760. 


The table presents the value of K for which g = 2Kp + 1 is the smallest divisor 
of M, rather than presenting the divisor, g. This has been done because, first, the 
value of K indicates more about the character of the divisor than g does, and 
second, to save space. All primes between 10,000 and 15,000 were examined. If 
any such prime is not listed in the table, it means that 2” — 1 has no prime factor 
<10-2”. 

Several criteria have been discovered concerning the divisors of the Mersenne 
numbers, M, . The best known of these (for K = 1) is due to Euler [2]. It states 
that if p = 4L + 3 and gq = 8L + 7 are both primes then g divides M, . For K = 3, 
Pellet was the first to state [3] that gq = 6p + 1 divides M, if g can be expressed in 
the form 4(2a + 1)’ + 27b’, and p and gq are both prime. For K = 4, Reuschle 
stated and Western proved [4] that g = 8p + 1 divides M, if ¢ can be expressed in 
the form d’ + 64(2c + 1) and p and gq are both prime. 

These calculations were performed on an IBM 650 system at Picatinny Arsenal. 
The program used was as follows: all prime factors g of M, (p > 2) are of the 
form q = 2Kp + 1, and of one of the two forms 8Z + 1. Thus, either K = 0 mod 
4 orp + K = 0 mod 4. Each prime p was expressed in binary form, p = >. ;-» a2", 
a; = Oor 1. The residues R; = Ri_, mod gq, Ro = 2, were found. Finally the residue 
[]i-o (R:)** mod g was evaluated. If this product is congruent to one (mod q) then 
q is a divisor of M,. 


Picatinny Arsenal 
Dover, New Jersey 


1... are ‘“Mersenne Numbers,’’ MTAC, v. 12, July 1958, p. 207 


2. G. Harpy & E. M. Wricut, An Introduction to the Theory of Numbers, 3rd edition, 
1956, Oniord University Press, p. 80. 
3. 


L. E. Dickson, History of the Theory of Numbers, v. I, Chelsea Publishing Co., 1952, p. 
Pd 


A. E. WEsTERN, ‘‘Some criteria for the residues of 8th and other powers, ‘“‘Proc. London 
Math. Sec. (2) 9 (1911) p. 244-272. 


On a Theorem of Mann on Latin Squares 


By R. T. Ostrowski and K. D. Van Duren 


Mann [1] proved the following theorem: Jf a Latin square L of order 4t + 2 
has a (2t + 1) X (2t + 1) block with as many as (2t + 1)* — t cells containing 
digits in a list of 2t + 1, then there exists no Latin square orthogonal to L. This theorem 
seemed until lately to give theoretical evidence for the truth of Euler’s conjecture 
that no pair of orthogonal Latin squares exists of any order 4¢ + 2. Now that Euler’s 
conjecture has been shown to be false [2], [3], a more detailed investigation of Latin 
squares of orders 4t + 2 seems worthwhile. 

The chief goal of the work reported in this note was to find an example indicating 
that Mann’s theorem is the best possible of its type for order 10—or conversely, to 
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accumulate experimental evidence suggesting that a more stringent theorem of the 
same nature remains to be proved. With the aid of the UNIVAC M-460 Computer, 
the former alternative was achieved. Displayed below is a pair of orthogonal Latin 
squares of order 10; the upper left 5 X 5 block of the first Latin square has 
5’ — 2 — 1 = 22 cells containing digits 0, 1, 2, 3, 4. (Note three italicized digits 
in each of the four separated 5 X 5 blocks of the left Latin square.) 


01234 56789 01923 84657 
34012 79865 67895 23104 
43120 97658 93746 58210 
12407 85396 38254 79061 
20375 68941 14507 36982 
57698 34120 25619 40873 
89756 12034 40138 62795 
65981 43207 56480 17329 
98563 01472 82071 95436 
76849 20513 79362 01548 


The second of the above Latin squares resembles either of the first pair of orthog- 
onal Latin squares of order 10 constructed by Parker [2]; it is disguised by permuta- 
tion of rows and columns. A UNIVAC M-460 program, coded by Parker [4], gener- 
ated the first of the above pair as an orthogonal mate of the second. Considering 
the results discussed near the end of this note it was largely luck that such an ex- 
ample as the above was found among the moderate number of cases examined. 

An observation shortened the computation time by a factor of four. In a Latin 
square of order 10 any 5 X 5 block contains each digit the same number of times as 
the block whose sets of rows and columns are complements of the sets of rows and 
columns for the first block. Further, the block with the same five rows and the com- 
plementary set of columns has the same sum of incidences for the set of five digits 
complementary to the set tallied in the initial block. Hence it is sufficient to inspect 
1(19)? — 15,876 blocks. This was implemented in the program by choosing all sets 
of five rows including the first, and all sets of five columns including the first. 

The authors’ program for the UNIVAC M-460 Computer carries out the follow- 
ing steps: 

1. Read into the computer a tape of one Latin square of order 10. 

2. Fora 5 X 5 block tally the incidence of each digit 0, 1, --- , 9. 

3. Sort out a set of five highest counts of the ten. Sum these five counts, and 
record this sum. If the sum for the block just run is the largest found for the Latin 
square, record the row and column indices. 

4. On completion of 15,876 (see above) passes through 2 and 3, generate an 
output tape. The information presented is the number of 5 X 5 blocks with sums 
of 25, 24, --- , 15; and the row and column indices for the first block giving the 
highest sum. 

5. Return to 1 or stop. 

Running time, including input and output, was about 40 seconds per Latin 
square. A repetitive process of this magnitude is barely thinkable without a com- 
puter; human scanning ability is not nearly good enough to locate a 5 X 5 block of 
highest sum of incidences of five digits. 
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The program was run for some sixty Latin squares of order 10 known to have 
orthogonal mates. Of these only one had a sum of 22—and this occurred for only 
one partitioning of rows and columns into fives. 

Another routine was written for the same computer to generate random Latin 
squares of order 10. Digits 0 through 9 were produced by an appropriately modified 
random number generator; cells were filled in sequentially subject to the required 
conditions, and changes made when completion became impossible. One hundred 
Latin squares were produced and processed by the first program. The highest sums 
of counts of five digitsin5 X 5 blocks were 19 for 76 Latin squares, 20 for 23 squares, 
and 21 for one; no Latin square of the hundred had a sum over 21. Of the hundred 
Latin squares, no two produced the same list of counts in 15,876 blocks; this evi- 
dence supports the claim of randomness. Before this study it had seemed plausible 
that most Latin squares of order 10 have sums over 22, since considerable searching 
by computer in recent years suggests that pairs of orthogonal Latin squares of 
order 10 are rare. Almost certainly there are necessary conditions for Latin squares 
of order 10 to possess orthogonal mates, aside from falsity of the hypothesis of 
Mann’s theorem. 


The problem was suggested by E. T. Parker. The authors thank him for help 
on various points. 


Control Data Corporation 
Minneapolis, Minnesota, and 
Remington Rand Univac Division 
Sperry Rand Corporation 

St. Paul, Minnesota 


1. Henry B. Mann, “On orthogonal Latin squares,” Bull. Amer. Math. Soc., v. 50, 1944, p 
249-257. 

2. E. T. Parker, “Orthogonal Latin squares,” Proc. Nat. Acad. Sci. U.S.A., v. 45, 1959, p. 
859-862. 

3. R. C. Boss, 8. S. SarrkHANDE & E. T. Parker, “Further results on the construction of 
mutually orthogonal Latin squares and the falsity of Euler’s conjecture,” Canad. J. Math., 
v. 12, 1960, p. 189-203. 


4. E. T. Parker, “A computer search for Latin squares orthogonal to Latin squares of 
order ten,” Notices Amer. Math. Soc., v. 6, 1959, abstract 564-71, p. 798. 


Complete Factorization of 2-1 


By K. R. Isemonger 


Algebraic factors of 2° — 1 are 2° — 1 = 7 and 2” — 1 = 6361-69431 -20394401, 
the latter factorization having been first published by F. Landry in 1869. 

R. M. Merson of Farnborough, Hants, England has found 13960201 to be a 
prime factor of (2"°° + 2° + 1)/7-6679. The resulting quotient is the integer 


N = 124 302 077 031 586 210 969, 


whose factorization was completed by me on 15 May 1960. 


The procedure employed was to exhibit N as the difference of two squares, 
namely, 





Received May 19, 1960; revised November 14, 1960. 








296 J. BOERSMA 


N=a -bD, 


where 


a = (2°-5-3°-53")a + 11 150 802 925. 
This representation of a can be deduced from theory presented by Kraitchik [1], 
combined with the fact that both —1 and 5 are quadratic residues of N, as estab- 
lished by suitable representations of N by quadratic forms. 
Corresponding to x = 102908, a” — N is the square of b = 114674787084. Hence, 


N is the difference of the squares of a = 115215488845 and of the preceding value 
of b. Thus 


N = 540 701 761-229 890 275 929. 
The primality of each of these factors was determined in a similar manner. The 
factorization of 2° — 1 is, therefore, now complete. 


19 Snape Street 
Kingsford, New South Wales 


1. M. Krartenrk, Théorie des Nombres, Gauthier-Villars et Cie, Paris, 1922, p. 146. 


Two Formulas Relating to Elliptic Integrals 
of the Third Kind 


By J. Boersma 


Using Legendre’s notation, the normal elliptic integral of the third kind is de- 
fined by the equation 


2 Pibes fetunkiy | ae tema pen : 
II(¢, ', &) -| (1 — a’ sin’ 0)+/1 — & sin? 6" 


For k’ < 1, the following expansion holds uniformly over the closed interval 


TT 
< sS-: 
0 0 5 





1 - -3 m7,2m + 2m 
/1 — FB sin? 6 = a ( Pe (—1) ke sin 6, 








ail ans A a NG dg cae a ak = 
where 7) = ber $2 Be orb 9. + Ot 8 form > 0, and *) =1. 
m m! 0 
The factor ae SS in the integrand is bounded for —2» <a < — 
1 — a? sin’? 0 sin? @ 


and 0 S @ S ¢; consequently, the expanded integrand may be integrated term by 
term. Such integration yields the series 


II (@, o’, k) = De bk, 
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where 


ies ee ease 
- m o l—a@sin?o ’ f 








° dé 1 ys 2 
hy eee eS ten ofr 
°= J i-eans Vig tan [1/1 — a tan ¢|,for — © <a <1, 


a 


= tan ¢, for a’ = 1, 





= : N4/e@ —1 2 l 
= Vani [Vo? — 1 tan ¢], for 1 < a < sas" 
In general, the coefficients b,, satisfy the recurrence relation 


2(m + 1)a’ays = (— 1)"*(2m + 1) (;) ton() + (2m + 1)Bu , 


m 
¢ 


where tem(@) = [ sin” 6 dé. 


0 
Byrd and Friedman [1] give [formula (902.00)] the recurrence relation 


2m — 1 ae 
tom(@) = ——— tom—2(p) — om sin” "¢ cos 


and explicit expressions for fo(@), fa(@), and &4(@). Corresponding to these we find 
bo —¢ 
20 


bs = ol [3a” sin ¢ cos @ + 6bo — 3(2 + a’ )¢] 
16a* 


by = 


5 


bs = [aR [2a' sin’ ¢ cos ¢ + a’ (3a’ + 4) sing cos + 8b — (8 + 40° + 3a‘)¢]. 


When ¢ = —x <a <1,andk < 1, we deduce the following expansion 


of the complete elliptic integral of the third kind: 


II(e’, k) =[] (5 a, t) = yi C, k™, 


bol 


where 











298 J. BOERSMA 
5x 2 4 8 
ae «san —_ 4 — —_ <~epeeammee . 
C3 | a 3a 8+ Vi | 


and, in general, the coefficients satisfy the recurrence formula 


2 
2(m + 1) emai = — (m + s) ~ Te + (2m + 1)em, 


which follows from the recurrence formula for b,, when use is made of the definite 


integral 
={% 1 1 
tom *) = [ sin’”6 dé = 1 ° T(m + 3)T() 
2 0 2 m! 


= 
_~ _1)"(-3 
=g(- 1) va . 

The expansions obtained above for Il (¢, a’, k) and [J (a’, k) constitute ex- 
tensions and simplifications of formulas (906.01) and (906.00), respectively, in the 
book already cited, by Byrd and Friedman. Furthermore, the coefficient of a” has 
been corrected here in the expression for c; appearing in (906.00). 


Mathematical Institute 
University of Groningen 
The Netherlands 


1. Paut F. Byrp & Morais D. Frrepman, Handbook of Elliptic Integrals for Engineers and 
Physicists, Springer-Verlag, Berlin, 1954. 
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REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 


62(C, D).—H. Scuuspert, R. Haussner & J. Ervesacnu, Vierstellige Tafeln und 
Gegentafeln fiir logarithmisches und trigonometrisches Rechnen in zwei Farben 
zusammengestellt, (Sammlung Géschen, Band 81) Walter de Gruyter & Co., 
Berlin, 1960, 156 p., 15 em. Price DM 3.60. 


This book, consisting largely of four-place tables, represents a revision by 
Joachim Erlebach of previous editions (dating back to 1913) by Robert Haussner 
of Schubert’s tables, which first appeared in 1898. 

As the full title states, the tables are printed in two colors—logarithms in red 
and all other material in black (in two previous editions that the reviewer has 
examined, the colors were brown and blue, respectively). A total of sixteen tables 
are given. These include 4D common logarithms of numbers from 1 to 10*, supple- 
mented by 7D logarithms of primes through 97; 4D common logarithms of trigono- 
metric functions; addition and subtraction logarithms; natural logarithms of inte- 
gers to 100, supplemented by 6D tables of the first ten multiples of In 10 and log e; 
squares and cubes; natural values of the trigonometric functions; values of fre- 
quently used mathematical constants, mainly to 6D; conversion tables (to 6D) 
from sexagesimal to radian measure; and miscellaneous useful tables, such as 
mortality tables, interest and annuity tables, geographical and astronomical tables, 
a table of physical properties of materials, and a periodic table of chemical elements. 

The table of antilogarithms appearing in earlier editions has been omitted; 
however, inverse interpolation in the table of logarithms is so simple, because of 
the small differences, that this deletion is of negligible practical importance. 

On the title page of the table of common logarithms appears a footnote explain- 
ing the formation of the logarithm tables of both numbers and trigonometric func- 
tions by truncation of the corresponding entries in one of the Bremiker editions of 
Vega’s tables. In cases of doubtful rounding to 4D of the Bremiker 7D logarithms, 
the convention was arbitrarily adopted of making the fourth-place digit odd. The 
reviewer found a total of ten entries in the table of logarithms of numbers that 
were subjected to this rounding; of these, four are affected by rounding errors. 
Such slight inaccuracies could have been avoided by referring to available ten- 
place tables such as those of J. Peters, which were first published in Germany in 
1922. These rounding errors appear also in the Schubert-Haussner tables dated 
1945; only one such error, however, appears in an edition dated 1917 that the re- 
viewer has examined. 

In summary, despite minor flaws the reviewer considers this book to be an un- 
usually extensive and valuable compilation of four-place tables. Their continued 
usefulness and popularity can be inferred from the numerous editions and printings 
of this member of the Sammlung Géschen series. 


J. W. W. 


63(C, L].—AnpbrEs Zavrotsky, “Construccion de una escala continua de las opera- 
ciones aritmeticas,”’ Revista Ciencia e Ingenieria de la Facultad de Ingenieria de 
la Universidad de Los Andes, Mérida, Venezuela, December 1960, No. 7, p. 38-53. 
The author studies the function S(z, y;n) = L~"(L"x + L"y), where Lz is 
the iterated logarithm: 
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L"x = log L” ‘x; L'x = log x. 


The function S(x, y; ) is of interest as a generalized arithmetic operation, 
since the values 0 and 1 of the parameter n yield x + y and «x-y, respectively. 

In a previous paper [1] tables for S(x, y; n) were given for xz, y = 0(1)10; = 
—1(1)3, where L”"x was defined in terms of logarithms to the base 2. 

In the current paper non-integer values of the parameter n are introduced by 
putting L"x = H(Gx — 1), where the mutually inverse operators H and G are 
defined by GLz = Gx — land H(x — 1) = LHx. 

In this paper, where Lx is defined in terms of natural logarithms, the function 
S(2, y;n) is tabulated for x, y = 0(1)10;n = 3, +/2, and for x, y = 2(1)10;n = -. 
All tabular entries are given to 5D. 


F. THEILHEIMER 
Applied Mathematics Laboratory 


David Taylor Model Basin 
Washington 7, D. C. 


1. A. Zavrotsxy, ‘‘Algunas generalizaciones del concepto de campo,’’ Acad. Ciencias Fis. 
Mat. y Nat., Caracas, Venezuela, Boletin, No. 28, 1946, 1947, 23 p. [See RMT 494, MT'AC, v. 
3, 1948-1949, p. 97.] 


64[F].—L. McKes, C. Nico. & J. Setrrinee, “Indices and power residues for all 
odd primes and powers less than 2000,’’ an unpublished mathematical table 
stored on magnetic tape, January 18, 1961. 


The computing center at the University of Oklahoma has recently computed a 
table of indices and power residues for all odd primes and powers thereof less than 
2000. The computations were done on a modified IBM 650 and have been stored 
on magnetic tape. Anyone desiring any portion of this table should write to: Direc- 
tor, Computing Center, University of Oklahoma, Norman, Oklahoma. 


AvuTHors’ SUMMARY 


65[G].—Ricuarp Bettman & MarsHauu Hatt, Jr., Editors, Proceedings of Sym- 
posia in Applied Mathematics, Vol. X, “Combinatorial Analysis,” American 
Mathematical Society, 1960, vi + 311 p., 26 cm. Price $7.70. 


This book contains the following papers, presented at a symposium on applied 
mathematics sponsored by the American Mathematical Society and the Office of 
Ordnance Research three years ago (April 1958). 


Marshall Hall, Jr. Current Studies on Combinatorial Designs 

R. H. Bruck Quadratic Extensions of Cyclic Planes 

D. R. Hughes On Homomorphisms of Projective Planes 

A. A. Albert Finite Division Algebras and Finite Planes 

L. J. Paige & The Size of the 10 x 10 Orthogonal Latin Square 
C. B. Tompkins Problem 

R. P. Dilworth Some Combinatorial Problems on Partially Ordered 

Sets 
R. J. Walker An Enumerative Technique for a Class of Combinatorial 


Problems 
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A. L. Whiteman The Cyclotomic Numbers of Order Ten 
A. J. Hoffman Some Recent Applications of the Theory of Linear In- 
equalities to Extremal Combinatorial Analysis 
A. W. Tucker A Combinatorial Equivalence of Matrices 
H. W. Kuhn Linear Inequalities and the Pauli Principle 
H. J. Ryser Compound and Induced Matrices in Combinatorial 
Analysis 
Marvin Marcus & Permanents of Doubly Stochastic Matrices 
Morris Newman 
A. M. Gleason A Search Problem in the n-cube 
D. H. Lehmer Teaching Combinatorial Tricks to a Computer 
J. D. Swift Isomorph Rejection in Exhaustive Search Techniques 
Olga Taussky & Some Discrete Variable Computations 
John Todd 
R. E. Gomory Solving Linear Programming Problems in Integers 
Richard Bellman Combinatorial Processes and Dynamic Programming 
Murray Gerstenhaber Solution of Large Scale Transportation Problems 
Robert Kalaba On Some Communication Network Problems 
J. D. Foulkes Directed Graphs and Assembly Schedules 
E. N. Gilbert A Problem in Binary Encoding 
M. M. Flood An Alternative Proof of a Theorem of Kénig as an 


Algorithm for the Hitchcock Distribution Problem 

Presumably the papers dealing with the application of computing equipment to 
attacks on combinatorial problems will be of most interest to readers of Mathematics 
of Computation. The editors list these as being the papers by Paige and Tompkins, 
Walker, Gerstenhaber, Flood, Gleason, Lehmer, Swift, Todd, and Gomory. Some 
of the papers listed are directly computational, that of Lehmer being an example, 
and some are more remote, Flood’s being not concerned directly with any computing 
instrument but highly influenced by his work on coding the assignment problem 
for a computer. 

The delay in publication has caused many of the papers to appear dated, but 
even these give some information not available in print elsewhere. The most spec- 
tacular example is the paper by Paige and Tompkins; this paper does not incorporate 
examples to be sought by methods of the paper, but obtained otherwise by work of 
several authors, R. C. Bose, 8. 8. Shirkhande, E. T. Parker, and others (mentioned 
in a footnote at the beginning of the paper). On the other hand, this paper outlines 
techniques which were then admittedly not feasibly applicable to the complete 
solution of the problem considered, but which are almost applicable on present day 
machines. 

Other papers, all of which were prepared by well chosen experts in their fields, 
seem to have greater durability, and this volume is one in which readers will find a 
great deal of current and useful material concerning combinatorial problems. On 
the whole, the volume is well worth perusing both for its computational implica- 
tions and for its general information concerning combinatorial problems. Marshall 
Hall’s opening paper, for example, presents an excellent survey of the field, which 
is valuable to almost anyone not completely up to date in combinatorial design. 
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The reviewer would like to note two improvements which would have greatly 
increased the value of the book. 

A unified bibliography at the end of the volume would be much handier than 
the separate bibliographies at the ends of the papers. Separate bibliographies are 
hard to find in the volume, and they are less useful for general reference to find a 
partially remembered work than a unified bibliography would be. Bibliographies 
associated with individual papers are somewhat more convenient to users of re- 
prints, but only the most insistent authors seemed to be able to talk the Society 
into furnishing reprints from this volume. On the whole, the reviewer believes that 
an integrated bibliography is almost a necessity in such integrated volumes, and he 
is continually puzzled when they do not appear. 

A more demanding task is the preparation of an instructive narrative summary 
of the type which has been produced in the Annals of Mathematics Studies devoted 
to games; [1] contains both a fine bibliography and a carefully written and highly 
lucid introduction. The editors of the present volume did make several interesting 
and instructive statements in their introduction, but they were not moved to shed 
the large amount of light they could have on the subject matter and its relevance. 

One interesting statement from the introduction might be worth quoting: 
‘What is very attractive about this field of research is that it combines both the 
most abstract and most nonquantitative parts of mathematics with the most 
arithmetic and numerical aspects. It shows very clearly that the discovery of a 
feasible solution of a particular problem may necessitate enormous theoretical 
advances. Perhaps the moral of the tale is that the division into pure and applied 
mathematics is certainly artificial and to the detriment of the enthusiasts on both 
sides. Furthermore, the way in which apparently simple problems require a complex 
medley of algebraic, geometric, analytic and numerical considerations shows that 
the traditional subdivisions of mathematics are themselves too rigidly labelled. 
There is one subject, mathematics, and one type of problem, a mathematical 
problem.” 

On the whole, the authors and the editors (the reviewer excepted) have done 
an excellent job of presenting the status of this rapidly developing field. The book 
is worth having if only for isolated papers: Lehmer’s, for example, showing how to 
make a computer behave, or the beguiling paper by Taussky and Todd, which 
attacks problems that would cause any computer and its masters to shudder. The 
editors, the American Mathematical Society, and the Office of Ordnance Research 
are to be congratulated on this volume. 

The reviewer, who usually shies away from reviewing his own work, undertook 
this review because of the difficulty of finding other non-contributing reviewers, 
and because of the dated and unimportant nature of the paper to which he con- 
tributed. 


C. B. Tompkins 
Institute for Defense Analyses 
Princeton, New Jersey 


1. A. W. Tucker & R. D. Lucs, Editors, Contributions io the Theory of Games, Volume IV, 
Ann. of Math. No. 40, Princeton University Press, Princeton, N. J., 1959. 
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66[G].—KennetH Horrman & Ray Kunze, Linear Algebra, Prentice-Hall, Inc., 
Englewood Cliffs, N. J., 1961, ix + 332 p., 24 em. Price $7.50. 


This book was written to provide a text for the undergraduate course in linear 
algebra at the Massachusetts Institute of Technology. It covers vector spaces, 
linear equations and transformations, polynomials, determinants, invariant direct- 
sum decompositions, the rational and Jordan canonical forms for matrices, inner 
product spaces, and bilinear forms. The treatment is imbued with the modern axio- 
matic and abstract spirit but many concrete illustrative examples and exercises 
are given so that a serious, persevering student can master it. When he has done 
so, he is well equipped to study not only applied mathematics and theoretical 
physics but also mathematical analysis. 

We heartily recommend the book either as a text or for private study. The 
authors are evidently teachers of experience and good judgment who know and 
like their subject. Their aim is high, and only the best students will master the 
course, but all the students will be the better for the part of it they learn. In these 
difficult days for students it is well to have texts so well planned and written as 
this one. 


F. D. MurNAGHAN 
Applied Mathematics Laboratory 


David Taylor Model Basin 
Washington 7, D. C. 


67(G].—G. Temp.e, Cartesian Tensors, John Wiley & Sons, Inc., New York, 1961, 
vii + 92 p., 19 em. Price $2.75. 


This introduction to vector and tensor algebra is well planned and should prove 
of value to the better than average student. The amount of material covered in less 
than 100 pages is surprising; in addition to the usual topics, there are chapters on 
isotropic tensors, spinors, and orthogonal curvilinear tensors. The influence of Weyl 
is evident in the treatment of isotropic tensors and of Brauer and Wey] in the treat- 
ment of spinors. We heartily recommend the book, which is addressed to first year 
students “pursuing an Honours course in Mathematics or Physics” in England. We 
found only two points where our instruction of the subject would differ slightly. 

First, in the reduction of a symmetric tensor to diagonal form the author dis- 


Saptatlg , . 
oon 8 and states that this must assume its lower bound 





cusses the ratio f(u) = —— 
when 0f/due = 0. This assumes, tacitly, that the lower bound is not assumed on 
the boundary of the region —1 S we S 1. The author’s statement that one must 
consider separately “the cases in which the matrices of S — \,u are of ranks 3, 2 
or 1” does not apply to Schur’s induction method, which is indifferent to possible 
equalities of the characteristic numbers of S. 

Second, in the treatment of spinors I would emphasize more the two-valued 
nature of the correspondence between rotations, R, and unimodular 2-dimensional 
unitary matrices, U. Thus, while to each U there corresponds only one R, to each 
R there correspond the two matrices + U. For me, the displacement, in modern 
physics, of vectors from their front rank position, in favor of spinors, is due to the 
fact that the 2-dimensional unimodular unitary group is not a proper representa- 
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tion, but rather an ambivalent one, of the 3-dimensional rotation group. Thus, 
there exist tensors of the unitary group which are not found amongst the tensors 
of the rotation group, and these tensors, the so-called spin-tensors of the rotation 
group, have their physical significance; on the other hand, all tensors of the rotation 
group are found among the tensors of the unitary group. 

F. D. MurNAGHAN 


68 [K].—E. D. Barractouan & E. 8. Paas, “Tables for Wald tests for the mean 
of a normal distribution,” Biometrika, v. 46, 1959, p. 169-177. 


Let x be normally distributed with unknown mean @ and known variance o’. 
The Wald sequential test for 6 = 6 against the alternative @ = 6,(@; > 9%) consists 
in taking observation 2,4; as long as a < (0, — 00)/ > fa: v;/o° + n(0> — 0,")/20° 
< b; sampling stopping when this relation first fails, with @ = 6 accepted if the 
left-hand inequality fails and @ = 6, accepted if the right-hand inequality fails. 
Let Z = —ac/(( — 0), h = (b — a)o/(0, — %), P_ = probability of accepting 
6 = 6) when true, P, = probability of accepting @ = @ when @ = 6,, N_ = average 
sample number when @ = 6), N. = average sample number when @ = @,, and 
uw = (0, — 6)/2. Table 1 of the Appendix contains 2D values of h and Z for P, = .05, 
-10(.10) .70, P_ = .95, .99, .995, .999, and w» = .25(.25)1.00. The values of a and b 
are determined by h and z. For Table 2 of the Appendix (the same combination of 
values for P, , P_, u occur as for Table 1), 2D values are given for N, and N_. 
Charts I-IV of the Appendix contain curves of P, and P_ as functions of h and Z 
for given u. These charts are obtained directly from Table 1 and can be used to 
determine the operating characteristics of the test for given a, b, o, 0, and @,. 
Also, a comparison is made between Wald’s approximations (to the operating char- 
acteristics and the average sample numbers) and the true values, for ten combina- 
tions of values for h, Z, u (Text-Table 1). The conclusion reached is that Wald’s 
approximations are not acceptably accurate for many applications when p 2 .25. 


J. E. Wats 
System Development Corporation 
Santa Monica, California 


69 [K].—A. T. Buarucna-Rem, Elements of the Theory of Markov Processes and 
their Applications, McGraw-Hill Book Co., Inc., New York, 1960, xi + 468 p., 
24 em. Price $11.50. 


In recent years the notions of probability have become increasingly important 
in the building of models of the world around us. This is true, for example, in certain 
of the physical sciences, social sciences, and in the simulation of military and other 
operations. Sometimes probabilistic notions appear directly as basic ingredients of 
the model, sometimes indirectly as the result of applying Monte Carlo methods to 
the solution of certain types of functional equations. 

The mathematical abstraction of an empirical process whose development is 
governed by probabilistic laws is known as a stochastic process. A special class of 
these processes are Markov processes in which the development subsequent to a 
time ¢t depends (probabilistically) only upon the state of the process at ¢ and not 
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upon its previous history. Their theory has been developed extensively in the past 
three decades and they have enjoyed increasingly wide application. Although the 
subject has been masterly treated from an advanced standpoint in Doob’s 
classical Stochastic Processes, there remains a dearth of textbooks in English which 
are accessible to the beginning graduate student who has not yet mastered the 
subtleties of modern analysis, especially measure theory. Notable exceptions, for 
the case of a finite number of possible states, are Feller’s highly regarded Introduc- 
tion to Probability Theory and its Applications and the recent book on Finite Markov 
Chains by Kemeny and Snell. 

The present work is intended as a graduate-level text and reference in applied 
probability theory. According to the author’s preface, its purpose is twofold: “first, 
to present a nonmeasure-theoretic introduction to Markov processes, and second, 
to give a formal treatment of mathematical models based on this theory which 
have been employed in various fields.”” Further, he states that the prerequisites are 
“a, knowledge of elementary probability theory (the first nine chapters of Feller, 
say), mathematical statistics (Mood level), and analysis (Rudin level). Some knowl- 
edge of matrices and differential equations is also required.” 

Part I develops the theory of Markov processes x(t), with particular reference 
to branching stochastic processes. Chapter 1, which lays the foundations for its 
successors, is devoted to the case wherein x assumes values in a discrete denumerable 
space, and the parameter ¢ is likewise discrete. In Chapter 2 the stochastic variable 
x again ranges over a discrete denumerable space, but the parameter ¢ is continuous. 
In Chapter 3 both z and ¢ are continuous. Each chapter is introduced by a summary 
paragraph. Thereafter the subject matter is developed carefully, systematically, and 
lucidly. Each chapter is followed by a list of problems and by a bibliography which 
is intended to be both complete and up-to-date. In order to make use of the latter, 
one would need a command not only of the standard languages of science—English, 
French, German, and Russian—but also of such less familiar languages as Polish 
and Hungarian. 

In Part II it is shown how the theory can be appled to problems in a variety of 
fields. These fields, together with a partial listing of the topics covered, are: 

Biology—Growth of populations, epidemics, theory of gene frequencies 

Physics—Cascade processes, radioactive transformations, particle counters, 
theory of nuclear reactors 

Astronomy and Astrophysics—Fluctuations of brightness in the Milky Way, 
spatial distribution of galaxies, radiative transfer 

Chemistry—Reaction kinetics 

Operations Research—Queueing theory and some of its applications. 

Extensive bibliographies are also given for each of these topics. 

The book closes with three appendices concerned with generating functions, 
Laplace and Mellin transforms, and Monte Carlo methods in the study of stochastic 
processes. 

The reviewer was very favorably impressed by this book. The quality of the 
printing is excellent, with pleasing page layout and very few printing errors. The 
author’s didactic skill and wide familiarity with the subject matter have resulted 
in a well motivated, well organized, clearly written text from which a student who 
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has the indicated prerequisites can readily gain an introduction to the theory and 
applications of Markov processes. 

R. P. Eppy 
Applied Mathematics Laboratory 
David Taylor Model Basin 
Washington 7, D. C. 





70 [K|.—R. E. Beckuorer, SALAH ELMAGHRABY & NorMAN Morsge, “A single- 
sample multiple-decision procedure for selecting the multinomial event which 
has the highest probability,” Ann. Math. Statist., v. 30, 1959, p. 102-119. 


Consider N k-nomial trials whose cell probabilities satisfy p, = --- = 
a. = pr./0 . We select that cell into which the most events fall, breaking a tie at 
random if it occurs. The authors give a 5D table of the probability of selecting cell 
k, for k = 2, 3, 4; 6* = 1.02(.02)1.1(.1)2(.2)3, 10; and N = 1(1)30. An approxi- 
mation is developed and compared with these values. 

J. L. Honess, Jr. 
University of California 
Berkeley, California 


71 (K].—K. G. Ciemans, ‘Confidence limits in the case of the geometric distribu- 
tion,’’ Biometrika, v. 46, 1959, p. 260-264. 


The author obtains confidence limits for estimating m, the expected number of 
trials before a device fails, given the sample mean £, and N, the number of devices. 
If N devices each are from an identical geometric distribution, the distribution of 
sample sums will follow a Pascal distribution. Two log-log charts are provided for 
two-sided 90%; and 98% confidence limits for m,1 S  S 10,000,and N = 2,5, 10, 
15, 20, 30, 50, 100. The charts are based on the exact distribution. For ¢ > 10,000, 
formulas and tables may be used to determine the confidence limits. For large 
N > 100 a special formula is given. Alternatively for large N, since sample means 
are approximately normal, confidence limits for m may be found as solutions of the 
quadratic equation obtained from t = +1/N(Z — m) + m(m + 1), where t is the 
usual normal deviate for the a percent point. 


L. A. AROIAN 
Space Technology Laboratories, Inc. 
Los Angeles, California 





72 [K|._E. T. Feperieut, “Extended tables of the percentage points of Student’s 
t-distribution,” J. Amer. Statist. Assn., v. 54, 1959, p. 683-688. 


The author states that in using Student’s ¢-distribution in testing component 
~ parts a need for extending the table of upper percentage points was revealed. The 
. method of calculation of these percentage points is presented, and a table containing 
these results is given. Let y, be the elementary density for Student’s ¢ with n degrees 


of freedom, and denote / y. dt by P. The values of t are given to 3D for P = 
to 


25, .10, .05, .025, .01, .005, .0025, .001, 5 X 10°, 25 X 10°, 1 X 10°, 5 X 10°, 
25 X 10°, 1X 10°, 5 X 10°, 25 X 10", 1 X 10°, 25 X 10°, 1 X 10°’, and n 
= 1(1) 30 (5) 60(10) 100, 200, 500, 10°, 2 x 10°, 10‘, and «. It would have been 
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advantageous had the large values of n been arranged conveniently for harmonic 
interpolation, such as n = 60, 120, 240, 480, 960, etc. 
L. A. AROIAN 


73(K].—_Irwin GuTrMan, “Optimum tolerance regions and power when sampling 
from some non-normal universes,” Ann. Math. Statist., v. 30, 1959, p. 926-938. 


This paper is concerned with obtaining 8-expectation tolerance regions which are 
minimax and most stringent (see [1] and [2]) for the upper tail of the single ex- 
ponential population and for the central part of the double exponential distribution. 
The single exponential probability density function (pdf) is of the form 
i exp [—(x — u)/o] with x = uy, where one or both of u and o are unknown. The 
double exponential pdf is of the form (20) exp (— | « — uw |/o), where uw is known 
and o is unknown. The sample values are 7 <--- < 2; € = Dime Ze n; 
$= ine (x; — a1)/(nm — 1); wo and oo represent known values of yu and a; 
t= > 71 | 2; — wo|. Then the optimum tolerance intervals, which are easily identi- 
fied with the situations considered, are [as(Z — wo), ©), [zi — Ddgoo, ~), 
[v1 — cgs, ©), and [wo — det, wo + det]. Tables I-IV contain 6D values of ag , bg , 
cs, dg , respectively, for n = 1(1)20, 40, 60 and 8 = .75, .90, .95, 99. The power 
of tolerance intervals is expressed in terms of parameter a, , where a; is determined 


as the solution of (ae) exp [—(x — u)/ac dx = y = measure of desirability, 
1(8) 


for the single exponential case, and from (2ae yf exp (—| x — uw | ac) dx = y for 
1(B) 


the double exponential case. Here /(8) is the tolerance interval considered and 
0 < y < 1 (large values indicate greatest desirability). Tables V, VI, and VIII 
contain 7D values of the power for intervals [ag(Z — wo), ©), [r1 — bgoo, ©), 
[uo — dgt, wo + dot], respectively, for n = 1(2)7, 10, 15, 30, 60, and 6 = .75, .90, 
.95, .99; likewise for x,cgs and Table VII, except that n = 2(2)10, 15, 30, 60. 


J. E. WatsH 


1. D. A. S. Fraser & Irwin Gutrman, ‘Tolerance regions,’’ Ann. Math. Statist., v. 
27, 1956, p. 162-179. 


2. Irwin Gutrman, “On the power of optimum tolerance regions when sampling from 
normal distributions,’’ Ann. Math Statist., v. 28, 1957, p. 773-778. 


74(K|.—Minos Jirexk & Orakar Likar, “Coefficients for the determination of one- 
sided tolerance limits of normal distribution,” Ann. Inst. Statist. Math. Tokyo 
v. 11, 1959, p. 45-48. 


It is well known that a random sample of size N from a normal universe with 
mean u and variance o° yields one-sided tolerance limits (— «, T,,) and (T,, +=) 
each of which includes at least a fraction a of the universe with probability P, where 


T. = £ + ke, 


T, =f-— ks 


, 
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and where, 


& 
ll 


n 
Zz zi/n, 
t=1 


[<) 


wn 
ll 


> (a: — €)*/(n — 1), 
i=1 
and 


/nk = t(n — 1, Uar/n, 1 — P). 


Here t(f, 6, €) [1] is the 100€ percentage point of the non-central ¢ distribution with 
f degrees of freedom, 6 is the measure of non-centrality in the definition of t, and u. 
is the 100(1 — a) percentage point of the unit normal distribution with zero mean. 

By use of the tables (especially Table IV) and iteration of the approximations 
given by Johnson and Welch in [1] the authors obtain values of the coefficient 
Vnk to 4S, for n = 5(1)20(5)50(10) 100(100) 300, for P anda = .90, .95, .99. 
A method for determination of these coefficients is given in [1], but the calculations 
are, of course, quite tedious, so that the present tables render a valuable service 
for practical applications to one-sided tolerance limits. 


S. B. Lirraver 
Columbia University 


New York, N. Y. 


1. N.S. Jonnson & B. L. Wetcu, “Applications of the non-central t-distribution,”’ 
Biometrika, v. 31, 1939, p. 362-389. 





75(K].—J. Pacuarss, “Tables of the upper 10% points of the Studentized range,” 
Biometrika, v. 46, 1959, p. 461-466. - 


Let q = w/s, where w is a sample range based on n values, and s is an inde- 
pendent estimate of standard deviation based on m values. Then tables of gq’ have 
been prepared for Pr(q 2 q’) = a, where a = 01, .05, .10, n = 2 (1) 20, and 
m = 1 (1) 20, 24, 30, 40, 60, 120, ©. Three significant figures are given throughout. 
The work of Harter [1] has been used in improving the accuracy throughout, par- 
ticularly for a = .01. For a = .01 and .05, these tables correct errors in [2]. 


I. R. SavaGe 
University of Minnesota 
Minneapolis, Minn. 


1. H. L. Harter, The Probability Integrals of the Range and of the Studentized Range, 
WADC Technical Report 58-484, vols. I & II, 1959. Office of Technical Services, U. S. Dept. of 
Commerce, Washington 25, D.C. 

2. J. M. May, “Extended and corrected tables of the upper percentage points of the 
‘Studentized’ range,’’ Biometrika, v. 39, 1952, p. 192-193. [RMT 1080, MTAC, v. 7, 1953, p. 94] 





76[(K].—K. C. S. Prmuar & Paso Samson, Jr., “On Hotelling’s generalization of 
T’,” Biometrika, v. 46, 1959, p. 160-168. 


Let S,/n; , S2/n2 denote independent covariance matrices arising from samples 

° ° : r(s) lo 
of sizes n, and n, from two p-variate normal populations, and U“’ = trace Sz S,, 
where s is the number of non-zero roots. Two approximations are compared with the 
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exact values for the upper 5 and 1 percentage points of U® for several values of 
= (m — 8 — 2)/2andn = (nm, — 8s — 2)/2. The approximations for the upper 5 

paps 1 percentage points of U” and U™ are given to 3 or 4D for m = 0, 5, 
= 15(5)50, 60(20) 100. 


I. OLKIn 
University of Minnesota 


Minneapolis, Minn. 


77(K).—J. G. Saw, “Estimation of the normal population parameters given a singly 
censored sample,” Biometrika, v. 46, 1959, p. 150-159. 


As estimators of the mean and variance of a normal distribution, given an 
ordered sample x, < 22 --- < 2, censored above z, , the author proposes 


r—l 
u* = Z1+ (1 —)z,, where Z,, = > ai/(r — 1), 


i=l 


r—1 


r—l1 
n* =a 2 (x; —2,)? +8 > (x; — 2)’, 


respectively, where ¢ is chosen to make u* unbiased, and @ and @ are chosen to make 
n* unbiased and of minimum variance. To facilitate use of these estimators, three 
tables are appended. Table 1 consists of entries of the weight factor « and Var (u*/c) 
to 10D for 1 < r < n S 20. Table 2 contains coefficients of (n + 1)~ in series 
approximations to « and to Var (u*/o). Weight factors a and 8 are not tabulated 
directly, and consequently routine application of the author’s estimates may be 
hampered. However, in order to permit calculation of these factors, Table 3, con- 
taining coefficients of (n + 1)~* in series approximations to them, has been in- 
cluded. These entries are given to 6D for p, = .50(.05).80, where p, = r/(n + 1). 


A. C. Conen, Jr. 
University of Georgia 
Athens, Georgia 


78[K].—Mrnorv Srorant, “The extreme value of the generalized distances of the 
individual points in the multivariate normal sample,” Ann. Inst. Statist. Math. 
Tokyo, v. 10, 1959, p. 183-208. 


Let xa’ = (Xia, *** »Xpa), @ = 1, --+ , nm, be m independent observations from 
a p-variate normal population with mean vector m’ = (m,, --- , m,) and covari- 
ance matrix A, and let @’ = (%, , --- , #,). The upper 5, 23, and 1 percentage points 
of the extreme deviate xmaxo = max,{(x; — #)’A'(x; — #)] is given to 2D for 


= 3(1)10(2)20(5)30, p = 2, 3, 4. When A is unknown, let L be a p X p matrix 
whose elements, /;; , are unbiased estimates of \;; , and have a Wishart distribution 
with v degrees of freedom. The upper 5, 23, and 1 percentage points of the Stu- 
dentized extreme deviate Tiaxo = max’ ((x; — #)'L~'(x; — )] is given to 2D for 
= 3(1)12, 14, » = 20(2)40(5) 60, 100, 150, 200. 
I. OLKIN 
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79[L].—L. Fox, Tables of Weber Parabolic Cylinder Functions and Other Functions 
for Large Arguments, National Physical Laboratory Mathematical Tables 
Volume 4, Her Majestv’s Stationery Office, London, 1960, iii + 40 p., 28 em. 
(Paperback) Price 12s. 6d. 


In a previous work [1], to tabulate, for instance, J,,(2) for large x, it was found 
convenient to write I,(x) = (24x) ~"e*F,,(x) and to tabulate the auxiliary function 
F,(«) for 1/x = z = 0(0.001)0.05. This device is very economical as compared with 
tabulation as a function of zx, and the F,,(2) entries are easily interpolated. Using 
the same idea, the present volume gives tables which supplement well-known tables 
of certain transcendentai functions as described below. An introduction describes the 
methods of computation. For each table second central or modified central differ- 
ences are provided. In some cases modified fourth-order central differences are also 
given. 

Table 1. The exponential integral 


z 


Ei(z) 


Il 
~ 
®& 
2 
~ 


= Bi —2) = / te dt 
Ei(x) = & F(z), z=n2. 


Thus, positive z relates to Ei(x) ; negative z, to Hi( —x). The function F is tabulated 
to 10D for +z = 0(0.001)0.100. 
Table 2. Sine and cosine integrals 


Si(z) = [ t" sin t dt, Ci(z) = | t cos t dt 
0 2 


Si(z) = 44 — Pcosx — Qsinz, Ci(xz) = Psinzx — Qecosz 


Values of P, Q are given to 10D for z = x ' = 0(0.001)0.100. 
Table 3. Airy integrals 
The notation follows Miller [2]. Let z = ir = jx - 


. —1/2 —1/4 — . —1/2_ —1/4 Y 
Ai(x) gg te et R, Bi(z) = w*?a “eS 

. —1/2 1/4 —§ oy 7” a eine 
A’i(x) eee Ww, B(x) = wre X 


Ai(—2x) + jBi(—2z) = 9? “*e-*(7P — Q) 


ll 
tol 


A'i(—2) + jB'i(—2) = Pre *@(GU — V), = jg = (-1)"?, 0 = E + 4. 


Ss 


The values of R, S, W, X, P, Q, U, V are tabulated to 10D for z = 0(0.001)0.050 
Table 4. The error integral 


i 2) z 
—4h72 = 2 -_ 2 
/ e*" dt = a4 'e*'S, | oe’ dt =x eT 
0 


The values of S and T to 10D are provided for x” = z = 0(0.001)0.010. Since 
error functions are often used in the form | Yaa dt, tables based on the latter 


representation should also have been prepared. 
Table 5. Factorial functions 
This is a table of the gamma function, its natural logarithm, and derivatives of 
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the latter. Let 
T(1 + 2) = (24)e%a**"f, 


In {T(1 + x)} = 3 ln (27) + (@ + 4) mz —2 +, 
ke 
F*’ (2) = ss {in T(1 + 2)}, F(z) = nz + fe, g'=z, 


For each table z = 0(0.01)0.10. Tabular value of f, xf2 , xF’, and 2°F” are given to 
10D; g, to 12D; 2°F’”, 2‘F™, to 9D. 

Table 6. Weber functions 

The notation follows Miller [3]. Let 


W(a,x) = (2k/x)""f cos x, W(a, —z) = (2/ka)’"f sin x 


s 


ad te W(a, —z) = —(x/2k)"’g sin y, 


qq W(a, 2) = — (kx/2)""g cos ¥, 


2 2 
x=eot+ie -alnt+intix, P=wt+ hr -—alnz+ he —ie 


z=n', k=(1+e")*—e", og. = Imm {T(} + ia)}. 


Values of f, ¢, g, w to 8D are tabulated for a = —10(1)10, z = 0(0.005)0.100. 
Values of k and ¢g are also provided. Table 6A gives 8D values of ~ for 
a = 0(0.05)2.50(0.1) 10.0. 


Y.L.L. 


1. L. Fox, A Short Table for Bessel Functions of Integer Order and Large Arguments, Royal 
Society Shorter Mathematical Tables, No. 3, Cambridge, 1954. See also Review 37, MT'AC, v. 
9, 1955, p. 73-74. 


2. J.C. P. Miuuer, The Airy Integral, giving Tables of Solutions of the Differential Equation 
y” = zy, British Association Mathematical Tables, Part-Volume B, Cambridge, 1946. See also 
Review 413, MTAC, v. 2, 1946-47, p. 302. 


3. NATIONAL Puysicat LaBoratory, Tables of Weber Parabolic Cylinder Functions. Com- 
puted by Scientific Computing Service Limited; Mathematical Introduction by J. C. P. Miller, 


Editor. Her Majesty’s Stationery Office, London, 1955. See also Review 101, MTAC, v. 10, 
1956, p. 245-246. 


80[X].—Germunp Dau tauist, “Stability and error bounds in the numerical in- 
tegration of ordinary differential equations,” Kungl. Tekn. Hégsk. Handl. 
Stockholm (Transactions of the Royal Institute of Technology, Stockholm, 
Sweden) Nr. 130, 1959, 85 p., 25 cm. Price Kr. 9. 
In the first part of his thesis, the author investigates stability of certain linear 
operators, 


(1) L = p(E) — hi’ (d'/dx’)o(E) + h'™ (d*"*/da™™)r(B), 


associated with numerical integration formulas for sets of differential equations of 
the form 


(2) d'g/dx’ = f(x, 9), 


where g and f are s-dimensional vectors, p(¢), o(¢), and 7(¢) are polynomials of 
degree k with real coefficients, and E is the displacement operator defined by 
Eu(x) = u(x + A), for any function u(x). The number k is called the order of L. 
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There exists an integer, p, called the degree of L, such that Lu(x) = O(h?*’) as 
h — 0. L is said to be stable if all zeros of p(¢) are of modulus <1, and the zeros 
with modulus 1 are of multiplicity Sr. 

The author shows that stability of L is equivalent to “stable convergence’’ of 
certain solutions of the difference equation associated with L to the solution of 
equation (2). If L is unstable the numerical solutions are “strongly unstable,” and 
the integration formula is practically useless. Several theorems are given concerning 
the largest possible degree of stable operators of a given order. For example, if 
7(¢) = 0, Lis unstable if p > 2({k/2] + [(r + 1)/2]); here [x] denotes the largest 
integer <2. 

If the initial-value problem, dy/dx = qy, y(0) = 1, ¢ constant, is treated nu- 
merically with a stable operator L of order k and degree k + 2, where r(¢) = 0, 
the solution is a linear combination of basic solutions, ¢,, 1 <j S k,n = 0,1, 
2, ---, where ¢;, are the roots of the characteristic polynomial, p(f) — qho(f). 
Let ¢; be the roots of p(¢); in the present case, | ¢;| = 1 and ¢; is single for all j. 
The author shows that, for h — 0, €j, ~ £;(1 + k,gh), where the k; = o(f£;)/ 
[¢;e’(¢;)] are called growth parameters and are real numbers, then 
tin ~ £;" exp (kjqhn). One of the £; , say {1 , is equal to 1, and k = 1. 

If Re(q) < 0, the solution ¢7, ~ e*” decays, but may be dominated by some 
of the other basic solutions which increase or decrease more slowly. If there exist 
solutions which increase for Re(q) < 0, the operator L is called weakly unstable; 
this occurs if at least one of the k; ,7 2 2, is <0. It is shown that a stable operator 
of even order k and maximum degree k + 2 is weakly unstable; such an operator 
generates an oscillatory solution whose amplitude increases at least as rapidly as 
| exp (—gnh/3) | if Re(q) < 0. 

In the second part of this thesis the author is concerned with estimating the 
norm of the error vector for r = 1. First, the error is évaluated for the linear vari- 
ational system d2/dx = B(x)zZ associated with (2), where B(x) is the Jacobian, 
(Of/8%)j—g2) - Then the effect of the linearization is estimated separately. Three 
error formulas are given. The second one, involving the directional derivative, 
p[B(x)] = limysoy [|| JZ + B(x) || —1], yields particularly good results if the 
numerical solution is smooth and if u[B(x)] < 0. 


WERNER LINIGER 
IBM Research Center 
Yorktown Heights, New York 


81[X]—A Koreanorr, with the collaboration of L. Bosserr, J. L. Gropomior & 
J. Jounson, Méthodes de Calcul Numérique, Tome I: Algébre non linéaire, Dunod, 
Paris, 1961, xxvii + 375 p., 25 em. Price 58 NF. 


The volume being reviewed is the first of a projected series. The characteristic- 
value problem is included, but matrix inversion is not, except to the extent that an 
iterative method that applies to a system of nonlinear equations would apply also 
to a system of linear equations. 

The French literature on numerical analysis is sparse indeed, and hitherto has 
been but slightly affected by the advent of the electronic computer. This book goes 
far toward filling the gap there, and would be a substantial contribution to the 
literature in any language. 
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A short initial chapter discusses digital computers and programming in rather 
general terms. Chapter 2, “Erreurs,” consisting of 70 pages, is perhaps the most 
complete and systematic single account in existence, discussing both rigorous bounds 
and probabilistic estimates. Chapter 3 considers iterative methods in general as 
applied to transcendental equations and to systems of equations. Chapter 4 is de- 
voted to polynomials specifically, and Chapter 5 to characteristic values. Each 
chapter is supplemented by an extensive bibliography. An appendix gives numerical 
results of applying the various methods to some specific problems. 

One can find fault with a few details. The method of Graeffe is given in the 
chapter on polynomial equations, although it is equally applicable to transcendental 
equations within a circle of analyticity. In the chapter on the characteristic-value 
problem, there is an interesting general section which derives the Jordan normal 
form, and develops the standard localization theorems of Gershgorin, Brauer, and 
others. But the treatment of the methods is somewhat disappointing, especially of 
the so-called “direct methods” (actually methods of reduction, although the re- 
viewer himself is guilty of having propagated the misnomer). It is even stated that 
the method of Krylov is not to be recommended for large matrices, whereas it has 
been shown by Bauer and this reviewer that most of the known direct methods are 
only specializations of the Krylov method. 

However, the authors could justify the rather sketchy treatment of the methods 
of reduction by arguing that these do not differ basically from methods of inversion, 
and are hence peripheral to their main interest, which would be the solution of the 
polynomial equation that can be obtained from a complete reduction. The treatment 
of this problem is by no means sketchy. In fact, it is the most complete and up-to- 
date account of known methods now in the literature, exhibiting well their inter- 
relations, and often presenting them from a fresh point of view. This volume can be 


recommended highly, and one can hope that the subsequent ones will be equally 
well done. 


A. 8. HousEHOLDER 
Oak Ridge National Laboratory 


Oak Ridge, Tennessee 


82[X|.—Hersert E. Sauzer, Dexter C. SHovuirz, & Exvizasern P. THompson, 
Tabies of Osculatory Integration Coefficients, Convair (Astronautics) Division of 
General Dynamics Corporation, San Diego, 1960, 43 p., 28 em. (Paperback) 


This publication contains tables to implement the use of quadrature formulas 
of the so-called osculating type, i.e., formulas in which the value of the integrand 


function and of its derivative at each point are used. The explicit formula con- 
sidered is 


rotph {n/2] 
(1) fi flu)du=h SOM (yf + DICH + Ral). 


o+ah i=—[(n—1) /2] 


The authors treatn = 2,3, 4,5;forn = 2 and 4, qis taken tobe 3, and forn = 3 
and 5, q is taken to be 0. These choices of g permit use of symmetry relations to 
reduce the amount of tabulation. The coefficients C{”(p) and DS” ( p) then are 


polynomials in p, and these are listed for each 7 indicated in (1). Ten-decimal tables 
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of C$” (p) and D$” (p) are given for the same values of i and for p = q(0.01){n/2], 
and nm = 2(1)5. Tables to aid in computation of an error bound are given. 

The authors state that these tables were calculated on an IBM 704 and listed 
to 12 decimal places. These were subjected to “functional” checks using a desk 
calculator, and then were rounded to 10 decimal places. On the basis of these checks 
the authors believe the coefficients are ail correct within about one unit in the tenth 
decimal place. 

It seems to the reviewer an unfortunate circumstance that the details of the 
checks made with the IBM 704 and the desk calculator were not more precisely 
given. It would have been simple, for example, to have checked the integrals of 
monomials in the range of precision, and this would have made an excellent inde- 
pendent check of accuracy, well worth the additional time required for automatic 
ealculation. 

The authors express great enthusiasm over the accuracy of the method and 
illustrate it with four examples which indicate rather well the situations in which 
osculating formulas may be used. These situations are ones in which the derivative 
of the integrand is available without excessive additional work. In particular, the 
suggested applications to orbit calculations in which the position and velocity 
vectors are known seemed very appropriate. However, this very example also sug- 
gests that there would be some real interest in extending the domain of p so that 
extrapolations could be made. In that case n = 1 would have been a possible choice. 

Mrs. Frieda Cohn of the Numerical Analysis Laboratory at the University of 
Wisconsin has calculated the integrals of a few monomials using these tables, and 
found that these checked within possibilities of round-off error. 

P. C. H. 


83[X, Z].—F. A. Ficken, The Simplex Method of Linear Programming, Holt, Rinehart 
& Winston, Inc., New York, 1961, vi + 58 p., 24 em. Price $1.50. 


The book literature on the subject of linear programming is in a state of rapid 
increase and Ficken’s contribution on the simplex method, the essence of the 
success of linear programming, is a welcome addition. It provides a rigorous treat- 
ment of old ideas in a 36-page discussion of duality, feasibility, boundedness, con- 
sistency, the simplex tableaux, degeneracy, etc. Brief material on inequalities, linear 
spaces, matrices, etc., appear in one appendix and theorems on existence and duality 
in a second appendix. The author’s bibliography makes no specific mention of G. 
Dantzig, to whom the topic of the book owes its existence and who has contributed 
immeasurably to the development of linear programming. A second point is the 
absence of material on integer programming. It seems important that a book on the 
subject, appearing in 1961, should shed light on this significant development. The 
author has avoided mention of this new gem in which the simplex process plays an 
important role. Otherwise, the book is recommended for the material it treats and 
for its clarity and rigor. 

Tuomas L. Saaty 


Office of Naval Research 
Navy Department 
Washington, D. C. 
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EpitorsaL Note: The author wishes to state that subsequent to the publication of this 
book he learned that the simplex method can be attributed to Dr. G. B. Dantzig. Accord- 


ingly, he wishes to extend appropriate acknowledgment along with sincere regret for this un- 
fortunate omission. 


84(Z].—C. J. Bouwxamp, A. J. W. Dutsvestiun & P. Mepema, Tables Relating 
to Simple Squared Rectangles of Orders Nine through Fifteen, Department of 
Mathematics and Mechanics, Technische Hogeschool, Eindhoven (Netherlands), 
August 1960, ii + 360 p., 19 cm. 


The problem of dividing a rectangle into a finite number of non-overlapping 
squares is associated with the problem of dividing a square into smaller squares, 
all different in size. The latter had been a fascinating but unsolved problem until two 
solutions were published independently. The first was by R. Sprague (Mathematische 
Zeitschrift, v. 45, 1939, p. 607) and the second by R. L. Brooks, C. A. B. Smith, 
A. H. Stone and W. T. Tutte (Duke Mathematical Journal, v. 7, 1940, p. 312-340). 
The latter paper accomplished this dissection into as few as 26 squares. Later papers 
by T. H. Willcocks (Fairy Chess Review, v. 7, 1948, and Canadian Journal of Mathe- 
matics, v. 3, 1951, p. 304-308) exhibited a dissection into 24 squares. It is not known 
whether a dissection into fewer squares exists or not. Part of the motivation for the 
present work is the search for such a dissection. 

Each of the foregoing dissections contained a subset of squares which formed a 
rectangle. A dissection of a rectangle which does not include a rectangular subset is 
said to be simple; and one in which no square is repeated is said to be perfect. 

Many dissections of rectangles have been derived and published in a series of 
papers. However, a systematic catalog of these dissections had been prepared only 
as far as 14 squares. C. J. Bouwkamp, who has been a principal investigator in this 
field, has collaborated with the other two authors to extend the work to 15 squares 
and to produce this book, which summarizes the known results. There are only two 
pages of explanatory text and the rest of the book is devoted to the tables. Table I 
is a list of the dissections of the 3663 simple perfect rectangles which contain up to 
15 squares. They are ordered by increasing values of the ratio of the sides of the 
rectangle. Table III contains the same rectangles, but they are listed by increasing 
values of the number of squares. Table IT lists the 342 pairs of rectangles which have 
the same ratio of sides. Of these, 13 pairs have no squares in common. Hence, a non- 
simple perfect squaring of a square is made by combining a pair of these rectangles 
with two squares whose edges are the dimensions of the rectangle. Table IV lists all 
of the 431 simple imperfect squared rectangles. Of these, four are square. 

Several peculiar properties of some of the items in these tables are mentioned. 
The computation of these rectangles was done on the IBM 650. The details of the 
analysis and the programming will be published later. 

One might consider the extension of these tables to 23 squares to settle the ques- 
tion of the smallest number of unequal squares which make a square. The reviewer 
estimates that this would require 16,000 books like this one! 


MIcHAEL GOLDBERG 
Bureau of Naval Weapons 


Washington, D. C. 
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85[Z|.—Food & Agriculture Organization of the United Nations, Handbook on 
Data Processing Methods, Columbia University Press, New York, 1959, vi + 111 
p., 27 cm. Price $1.00 (Part 1, Provisional Edition) Paperback. 


This handbook is a well-written elementary discussion of the methods of process- 
ing statistical data from surveys or censuses. It reviews the various methods of data 
processing, with emphasis on hand and punch card methods. Only one page is 
devoted to stored program computers. A bibliography is included with the discussion 
of each topic. The general principles concerned with planning, organizing, adminis- 
tering, and operating data-processing services are adequately covered. Detailed 
instructions are given for the proper handling and preservation of records and for 
sorting records, using manual or punch card methods. The handbook is intended for 
statisticians or supervisors of data processing activities, and hence does not cover 
the special operating features and details of available machines. Although prepared 
primarily for the use of statisticians in countries unfamiliar with present data- 
processing tools and procedures, it could serve as a useful introduction to the subject 
for others. 

A. EUGENE SMITH 
Bureau of Ships 


Navy Department 
Washington 25, D. C. 


86(Z]—Fritz Rupotr Gintscu, Hinfiihrung in die Programmierung Digitaler 
Rechenautomaten, Walter de Gruyter & Co., Berlin, 1960, 144 p., 24 cm. Price 
DM 24. 


This introduction to the programming of digital computers assumes very little 
previous knowledge of computational procedures. The various steps in the prepara- 
tion of a computer program are explained with one particular computer, the Z22, 
in mind. A large number of specific problems are worked out in considerable detail 
by giving the algorithm, the flow chart, and the detailed coding. Much attention is 
paid to the discussion of loops, and subroutines are treated with considerable detail. 
Since the machine code of the Z22 contains only very few instructions—not even a 
multiplication instruction—nearly all practical coding has to be done in an external 
code. Therefore, all the illustrative coding in the book up to the second-to-the-last 
chapter is done in an external code. 

The book is devoted exclusively to introducing the reader to programming. It, 
therefore, does not contain any chapters on numerical analysis or the details of 
machine logic. 

F. THEILHEIMER 


87[Z].—-ALLEN Kent, Editor, Information Retrieval and Machine Translation, Part 
1, Interscience Publishers, Inc., New York, 1960, xv + 686p., 23 em. Price $23.00. 


This book contains twenty-one of the papers presented at the International 
Conference for Standards on a Common Language for Machine Searching and 
Translation, sponsored by Western Reserve University and the Rand Develop- 
ment Corporation, and held in Cleveland, Ohio, September 6-12, 1959. It is the 
first part of Volume III (in two parts) of the series entitled Advances in Documenta- 
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tion and Library Science. The papers, plus discussions and numerous references, 
afford a good indication of work now in progress on mechanizing the storage and 
retrieval of information, particularly as it relates to scientific documentation. 

The first article (by Allen Kent) reviews the current state of the art in machine 
literature-searching and translation. It includes in tabular form a comprehensive 
survey of searching devices and analytical techniques under investigation. Addi- 
tional tables provide a significant sampling of the recent literature bearing on 
machine translation. These are arranged according to institution and investigator, 
subject field, languages involved, and type of equipment used. It is interesting to 
note that either English or Russian is the target language in nearly every one of the 
more than one hundred different investigations cited, while the most frequent 
source languages are Russian, English, French, German, Chinese, and Japanese, 
in approximately that order of frequency. 

The remainder of the book is devoted to topics connected with information 
retrieval. The papers are by workers from such fields as library science, patent law, 
engineering, chemistry, medicine, the behavioral sciences, linguistics, symbolic 
logic, and computer technology. The majority of them deal with problems in 
organizing, indexing, and coding information for mechanical processing. 

The proclivities of each discipline are revealed by the varied approaches to the 
problems. The librarian is interested in enumeration of characteristics which serve 
to describe and classify things; the patent searcher is concerned primarily with 
structure and function, and he needs elaborate means for specifying interrelation- 
ships among assemblages of components; the chemist must be able to search in 
terms of properties, composition, and arrangement; the medical researcher struggles 
to elucidate information through techniques of statistical analysis and correlation 
of recorded observations; the linguist has to do with syntactic and semantic analysis 
of messages to discover the intentional as well as intensional implications which they 
convey; the logician manipulates information in symbolic form and tries to form- 
ulate a consistent set of rules for drawing valid inferences and selecting optimum 
strategies in a given frame of reference; and the systems designer seeks to provide 
a means for satisfying the multitudinous requirements of the potential users on an 
economic basis. Such a diversity of emphasis is characteristic of a field which is 
evolving rapidly. 

Nevertheless, the problems of indexing and searching technical literature, 
patents, chemical structural formulas, medical records, etc., do have a common 
basis, even though they may differ in a number of nonessential features. There is 
general agreement that some kind of coordination of descriptors is necessary for 
efficient retrieval, but that something more sophisticated than simple conjunction 
of characteristics will be required to cope with the increasing size and complexity 
of subject matter files. More particularly, it has been found in encoding subject 
headings that the relations expressed by the product, sum, and difference of classical 
logic are not sufficient by themselves for achieving the desired degree of specificity. 
A reasonably complete representation also requires a means for denoting the syntac- 
tic relations between descriptors. 

A number of schemes are discussed, particular importance being attached to 
the rules for analyzing subject matter and to the notation for expressing the results 
of the analysis. These include the use of modulants, interfixes, and role indicators to 
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specify the types of relationships between terms and the use of phrase boundaries, 
parentheses, or other special marks of punctuation to segregate groups of asso- 
ciated terms. Such devices greatly increase the discriminating power of an index by 
preventing false associations of terms. 

However, there is an evident lack of agreement on the means for attaining these 
desirable ends. This stems from difficulties inherent in linguistic analysis. It is neces- 
sary to discover and to define the syntactic relationships which are significant in 
mechanical indexing and searching. These are sometimes elusive, although such 
categories as species, genus, part, whole, factor, process, product, and attribute are 
widely recognized. The problem of synonyms (or near synonyms) and recognition 
of specific-generic connections can best be handled through the preparation of 
technical dictionaries and thesauri. These will provide a semantic analysis of com- 
plex notions in terms of more elementary concepts. 

Some concern is voiced that the potential gain in accuracy of representation 
may be offset by an increase in ambiguity of indexing as a result of the greater com- 
plexity of the task. No procedure can be considered satisfactory unless there is a 
reasonable prospect that two different workers indexing an item will arrive at sub- 
stantially the same description. Further research on the development of informa- 
tion retrieval languages is needed. 

The remaining articles describe the capabilities and limitations of several existing 
or proposed systems. They discuss the types of information services provided, the 
way in which the files are set up, and the types of mechanical devices used. Equip- 
ment mentioned includes Minicards, Magnacards, the WRU Searching Selector, 
the IBM 7090, and the new GE-225 computer system with special features for per- 
forming searches. 

The book is recommended to those who are currently working in some area 
touching upon information retrieval and to any others with a special interest in 
this fast-growing field. 

Tuomas 8. WALTON 
Applied Mathematics Laboratory 
David Taylor Model Basin 
Washington 7, D. C. 
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TABLE ERRATA 


301.—A. A. Bennett, W. E. Minne & H. Bareman, Numerical Integration of 
Differential Equations, National Research Council, Bulletin No. 92, Washington, 
D. C., 1933. Reprinted by Dover Publications, Inc., New York, 1956. 


On Page 83, in formula (14) the last term should read —yz},V°u, instead of 
tzipV'Un e 
H. E. Sauzer 
Convair-Astronautics 
San Diego, California 


302.—R. A. Buckincuam, Numerical Methods, Pitman Publishing Corporation, 

New York, 1957. 

On page 153, the following errors are noted. In Table 6.1 in the expansion of 
(2/6) sinh'(6/2)’, corresponding to r = —4 the numerator of the coefficient of 
5° should read —8 instead of —7. In the third line from the bottom of the same page, 
in the formula for &°J’u(1/2), the coefficient of 5° should read —367/193536 instead 
of —367/193537. In the last line, in the formula for &*I*u(1/2), the coefficient of 5 
should read 7/5760 instead of 5/5760. 

H. E. Sauzer 


303.—R. Ecrersp6rrer & L. Ecerspérrer, Formeln und Tabellen der zugeordneten 
Kugelfunktionen I. Art von n = 1 bis n = 20. (Reichsamt fur Wetterdienst, 
Wiss. Abh., 1(6).) Springer-Verlag, Berlin, 1936. 


The values R,,’ have been independently computed to 88 for all integers j, n 
such that n + 7 is odd and 0 S j < n S 19. Comparison of these values with 
corresponding entries in the above tables revealed just one error; namely, in the 
fifth line of the value of Rio, read +0.0067029 2050760 sin 10¢ instead of 
+0.0067028 2050760 sin 108. 

Frank K. BAMBERGER 
Operations Analysis Laboratory 
The University of Chicago 
Chicago 37, Illinois 


304—A. Erpféty1, W. Maenus, F. Operuertincer, & F. Tricomi, Tables of 
Integral Transforms, Volume 2, McGraw-Hill Book Co., Inc., New York, 1954. 


On page 49 the following errors have been noted. 
Equation (11) is in error unless n = 0. The right-hand side should read 
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Equation (13) is in error for negative values of n. In this case the right-hand 
side should read 


(-1)"8" "ry ae .(a8)K,(yB) 
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Equation (15) is in error for negative values of the integer n. In this case the 
right-hand side should read 


aie... 





42 n+1) — —(m+n+1) — 
= re mae SD apes (0) 
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Equation (16) is in error for negative values of n when 0 < y S a and for 
positive values of n when a S y < o. For these cases the right-hand side should 
read 
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The following errors have been noted in section 16.2, pages 277-279. Equation 
(10) should read 


[ (z — x)*2"P, (x) dx = 22"Q,(z); 
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equation (18) should read 


[ (2 — z)P,(2)P,(x) dz = 2P,(z)Qn(2); 


equation (25) should read 
1 
[ (z — 2) "(1 — 2°)””’P,"(x) dx = 2(2 — 1)”"Q,"(z); 


and equation (26) should read 


1 


[ a’ (z—2)"*(1 — a*)”"P,"(2) de = 22(2 — 1)”"Q,"(z). 


J. N. Newman 
Hydromechanics Laboratory 


David Taylor Modei Basin 
Washington, District of Columbia 


Eprtrortau Norte: Professor Erdélyi reports that the errors noted on page 49 were commu- 
nicated to him in 1957 by J. W. Stuart. 


J. W. W. 


305.—H. Scnusert, R. Haussner & J. Ertesacn, Vierstellige Tafeln und Gegen- 
tafeln fiir logarithmisches und trigonometrisches Rechnen ..., Walter de Gruyter 
& Co., Berlin, 1960. [See RMT 62, Math. Comp., v. 15, 1961, p. 299] 


The following rounding errors were discovered in Table I, “Vierstellige deka 
dische Logarithmen der Zahlen von 1 bis 10 000”: 


for read 
log 2298 3613 3614 
log 3560 5515 5514 
log 4023 6045 6046 
log 4719 6739 6738 


d. We We 
306.—T. N. Turexe, Jnterpolationsrechnung, Teubner, Leipzig, 1909. 


On page 96, in the formula for D~™, the coefficient of A°Vz’ should read 15/5760 
instead of 15/5790. 

On page 97, in the formula for D™, the coefficient of a°O.V,. should read 24/5760 
instead of 24/5770. 


H. E. Sauzer 
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